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We describe the current status of our numerical simulations for the collapse of a massive 
stellar core to a black hole (BH) and the merger of binary neutron stars (BNS), performed 
in the framework of full general relativity incorporating finite-temperature equations of state 
(EOS) and neutrino cooling. 

For the stellar core collapse simulation, we present the latest numerical results. We 
employed a purely nucleonic EOS derived by Shen et al. (Nucl. Phys. A 637 (1998), 435). 
As an initial condition, we adopted a 100 Mq presupernova model calculated by Umeda and 
Nomoto (Astrophys. J. 637 (2008), 1014), which has a massive core (M « 3M©) with a 
high value of entropy per baryon (s ~ 4fcs). Changing the degree of rotation for the initial 
condition, we clarify the strong dependence of the outcome of the collapse on this. When 
the rotation is rapid enough, the shock wave formed at the core bounce is deformed to be a 
torus-like shape. Then, the infalling matter is accumulated in the central region due to the 
oblique shock at the torus surface, hitting the proto-neutron star and dissipating the kinetic 
energy there. As a result, outflows can be launched. The proto-neutron eventually collapses 
to a BH and an accretion torus is formed around it. We also found that the evolution of the 
BH and torus depends strongly on the rotation initially given. 

In the BNS merger simulations, we employ an EOS incorporating a degree of freedom for 
hyperons derived by Shen et al. (Astrophys. J. Suppl. 197 (2011), 20), in addition to the 
purely nucleonic EOS. The numerical simulations show that for the purely nucleonic EOS, 
a hypermassive neutron star (HMNS) with a long lifetime (S> 10 ms) is the outcome for the 
total mass M < 3.0Mq. By contrast, the formed HMNS collapses to a BH in a shorter time 
scale with the hyperonic EOS for M > 2.7 Mq. It is shown that the typical total neutrino 
luminosity of the HMNS is ~ 3-10 x 10 53 ergs/s and the effective amplitude of gravitational 
waves from the HMNS is 2-6 x 10~ 22 at / » 2-2.5 kHz for a source distance of 100 Mpc. 

§1. Introduction 

Along with the development of formulations and numerical techniques, as well 
as progress in computational resources, numerical relativity (NR) is now the most 
viable approach for exploring phenomena accompanying strong gravitational fields, 
such as gravitational collapse of massive stellar core to a black hole (BH) or a neutron 
star (NS) and coalescence of compact-star binaries. These phenomena show a wide 
variety of observable signatures, including electromagnetic radiation, neutrinos, and 
gravitational radiation, and observations of neutrinos and gravitational radiation will 
provide us unique information of strong gravity and properties of dense nuclear mat- 
ter otherwise cannot be obtained. Next-generation kilo- meter-size gravitational- wave 
detectors such as LIGO, 1 ' VIRGO, 2 ' and KAGRA® will report the first detection 
of gravitational waves in the next ~ 5 years. In addition, the above phenomena are 
promising candidates of the central engine of long gamma-ray bursts (LGRB) and 
short gamma-ray bursts (SGRB).® 
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All four known forces of nature are involved and play important roles in the 
stellar core collapse and the merger of binary compact objects: General relativistic 
gravity plays a crucial role in the formation of a BH and a neutron star. Neutrinos 
produced by weak-interaction processes govern the energy and chemical evolution 
of the system. The electromagnetic and strong interactions determine the thermo- 
dynamical properties, in particular equation of state (EOS) of the dense nuclear 
matter. Strong magnetic field, if it is present, can modify the dynamics of the mat- 
ter motion. To study the dynamical phenomena in general relativity, therefore, a 
multi-dimensional simulation incorporating a wide variety of physics is necessary. 

We performed a simulation of stellar core collapse to a neutron star^ and a 
black holep2J incorporating a finite-temperature, a self-consistent treatment of the 
electron capture, and neutrino cooling by a detailed leakage scheme, for the first time. 
Such multi-dimensional simulations had not been done in full general relativity until 
quite recently Rl Ott et alP (see also Dimmelmeier et alP) performed fully general 
relativistic simulations of stellar core collapse, employing a finite-temperature EOS 
derived by Shen et alP 1 (Shen-EOS) for the first time. In their calculation, however, 
the electron capture rate was not calculated in a self-consistent manner and neutrino 
cooling is not taken into account. Instead, they adopted a simple parameterized 
prescription proposed by Ref. 1 10p : The electron fraction is assumed to be a function 
of density which is presumed based on a result of a single core collapse simulation 
with a specific initial condition. Recently, Miiller et al.^ performed simulations of 
stellar core collapse with detailed microphysics and neutrino transfer. However it is 
done in the framework of an approximate general relativistic gravity.^ Kuroda et 
alPS 1 have made a fully general relativistic code with an approximate treatment of 
neutrino transfer, applying the schemes developed by the authors-H^'EO 1 Ott et alP^ 1 
performed simulations of rotating stellar core collapse, employing the parameterized 
prescription^ 1 in the collapse phase, and a ray-by-ray neutrino leakage scheme after 
bounce. 

As for the compact-star binary mergers, there have been only a few studies in 



general relativistic frameworks In the framework of an approximate general 
relativistic gravity (the conformal flatness approximation^ 1 ) , Oechslin and Janka-^ 1 
performed simulations of binary neutron star (BNS) mergers adopting the Shen-EOS 
and a finite-temperature EOS by Lattimer and Swestyp^ 1 but they did not take 
account of weak interaction processes. Duez et alP^ 1 studied effects of EOS on the 
dynamics of black hole-neutron star mergers (BHNS) adopting the Shen-EOS in full 
general relativity. Recently, Bauswein et alM^* performed BNS simulations adopting 
a wide variety of EOS in the conformal flatness approximation and investigated the 
dependence of gravitational wave spectra on EOS. In these works, however, the weak 
interaction processes are not included. 

In this paper, we describe our latest results of numerical-relativity simulations 



*' There are a number of simulations of stellar core collapse in spherical symmetry^ in which 
Boltzmann's equation is solved and detailed microphysical processes are implemented. 

**' There are several studies of binary neutron star mergers in Newtonian frameworks in which 
finite-temperature EOS and weak interactions arc taken into account together with neutrino cool- 
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for the stellar core collapse to a BH 2< 9 and the BNS merger VdJtesJ wri ich are per- 
formed incorporating both a finite-temperature EOSPJ'^ and neutrino coolingP^ 
For reviews on other topics, namely, simulations of stellar core collapse to a neu- 
tron star, of BHNS merger, and of BH-BH binary merger, the reader may refer to 
Refe-BDH ,12211 JE}, Refs.[32D,[33D and Refs. EH) ,[35}, respectively. 

The paper is organized as follows. In § [21 we first briefly summarize basic 
equations, input microphysics, and numerical setup. The results of simulations of 
the stellar core collapse and BNS merger are described in § [3] and § [U respectively. 
Section[5]is devoted to a summary. Throughout this paper, H, ks, c, and G denote the 
Planck's constant, the Boltzmann's constant, the speed of light, and the gravitational 
constant, respectively. In appendices, details of the microphysics adopted in our 
latest implementation are summarized for the purpose of completeness. We adopt 
the geometrical unit c = G = 1 in § 12.11 and § 12.21 

§2. Basic Equations and Numerical Method 

2.1. Einstein's equations and gauge conditions 

The standard variables in the 3+1 decomposition of Einstein's equations are the 
three-dimensional metric 7jj and the extrinsic curvature Kij on a three-dimensional 
hypersurface defined by^ 

ItiP = 9^u + n^n u , (2-1) 
1 

where is the spacetime metric, is the unit normal to the three-dimensional 
hypersurface, and ££ n is the Lie derivative with respect to the unit normal . Then 
the line element is written in the form 

ds 2 = -a 2 dt 2 + j ij {dx i + f3 i dt){dx j + P j dt), (2-3) 

where a and fi % are the lapse function and the shift vector, which describe the gauge 
degree of freedom. 

Numerical simulations are performed in the so-called BSSN-puncture formula- 
tionpS | '!23 | 'H2l | i n which the spatial metric, is conformally decomposed as 7jj = 
W~ 2x fij where the condition, det(7jj) = 1, is imposed for the conformal spatial met- 
ric 7™. From this condition, the conformal factor is written as W~ e = det(7jj). The 
extrinsic curvature, K^, is decomposed into the trace part, K, and the traceless 
part, Aij, as = A{j + (l/3) r jijK. The traceless part is conformally decomposed 
as = W~ 2 Aij. To summarize, the fundamental quantities for the evolution 
equation are now split into W, K, and Ay. Furthermore, the auxiliary variable 
F{ = S^d^ij is introduced in the original version of the BSSN formulation. 3 ^ Mer- 
its of using W as a conformal factor are that (i) the equation for the Ricci tensor is 
slightly simplified, (ii) no singular term appears in the evolution equations even for 
W — > 0, and (hi) the determinant of jij is always positivePD'^ 
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The basic equations to be solved are 

d t - P k d k ) W=^(aK- d k (3 k ) W, (2-4) 

2 

F- 

1 



fa - P k 8 k ) % = -2aA i3 + j ik d 3 p k + j jk diP k - ^ijd k p k , (2-5) 



d t -P k d k )K = -D k D k a + a 



AijA lJ + -K 2 



+ 4to ( e£ otal + S Total ) , (2-6) 



8 t - (3 k d k ) A i3 = aW 2 (r i3 - \%R \ ~ (w 2 DiDja - ^ ij D k D k c^j 

+a (KA i3 - 2A ik A k ^ + A ik d 3 (3 k + A jk d^ k - ^A i3 d k p k 

-Sira (w*S(§** - l%S Tot ^ , (2-7) 

(dt - P k d k ) Ft = -16najJ otal 

+2a U kj d 3 A ik + Aikdjf* - \^ l dih 3l - 3A\d k In W - -fyA 

+V k {-2A l3 8 k a + (dkP^dihij 

+d k (ludjP 1 + JjAP 1 - pub/A | , (2-8) 

where = ^ — 5^ . ^R, ( 3 >Ri 3 , and Di are the Ricci scalar, the Ricci tensor, and 
the covariant derivative associated with three-dimensional metric 7«, respectively. 
The matter source terms are the projections of the stress-energy tensor (see Eq. 
(JESS)) with respect to and 7/u/ , and S™ = j»S^ taX : 

e Totai^ (T TotSLl ) a,3 n a ni3, (2-9) 

^-Total = _ (r Total )Q/ 9 7 . Qn/3) (2 . 1Q) 
^Total- (r Total )Q /9 7 . Q7 .^ (2 . n) 

where (T^° ) a R is the total energy-momentum tensor (see Eq. (|2T4p for definition). 
As a gauge condition for the lapse, we use a dynamical slicing^'® 

(d t - P k d k )a = -2Ka. (2-12) 

It is known that this dynamical slicing enables to perform a long-term-evolution 
simulation of neutron stars and BH spacetime. The shift vector is determined by 
solving a dynamical-shift equation^* 

d t p k = t\F l + Atd t F l ). (2-13) 

Here the second term in the right-hand side is necessary for the numerical stability, 
and At denotes the numerical timestep. 



*^ In the axisymmetric stellar core collapse simulation described in § [3] we do not include the 
advection term f3 k dka. 
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A fourth-order-accurate finite differencing in space and a fourth-order Runge- 
Kutta time integration are used in solving Einstein's equations and the gauge con- 
ditions. 

2.2. Hydrodynamic equations and GR leakage scheme 

Numerical simulations were performed using a fully general relativistic hydrody- 
namic code^-0 recently developed, in which a nuclear-theory-based finite-temperature 
EOS, a self-consistent treatment of electron and positron captures, and neutrino 
cooling by a general relativistic leakage scheme, are implemented. 

Because the characteristic timescale of the weak-interaction processes (twp ^ 
\Y e /Y e \) is much shorter than the dynamical timescale, idyn> in hot dense matters, 
source terms in the hydrodynamic equations become too stiff for the equations to 
be solved explicitly in a straightforward mannerPS A very short timestep (At < 
t wp <C tdyn) will be required to solve the equations explicitly. The characteristic 
timescale, ti ea k> with which neutrinos leak out from the system, by contrast, is much 
longer than t wp in the hot dense matter region, as ii ea k ~ L/c ~ tdyn, where L is the 
characteristic length scale of the system. Using this fact, we developed a method of 
solving the hydrodynamic equations in which the source terms are characterized by 
the leakage timescale £i ea k- 

Note that neutrino heating is not included in the current version of the leakage 
scheme. A conservative shock capturing scheme^ 1 with third-order accuracy in space 
and fourth-order accuracy in time is employed for solving hydrodynamic equations. 
In this section, we adopt the geometrical unit c = G = 1. 

2.2.1. Energy- momentum conservation equation 

The basic equations of general relativistic hydrodynamics including the radiation 
transfer for neutrinos are 



where (r Total ) a) g is the total energy-momentum tensor, and (T F ) Q p and {T L ') a p are 
the energy-momentum tensor of fluids and neutrinos, respectively. Equation (|2-14|) 
can be decomposed, by introducing the interaction source term, as 



Here the source term Q a is characterized by t wp and becomes too stiff in hot dense 
matter regions. To overcome the situation, the following procedures are adopted. 

1. The neutrino energy-momentum tensor is decomposed into 'trapped- neutrino' 
((T u > T ) a p) and 'streaming-neutrino' ((T u,s ) a /s) parts as 




(2-14) 




(2-15) 
(2-16) 




(2-17) 



Here, the trapped-neutrino part phenomenologically represents neutrinos which 
interact sufficiently frequently with matter, and the streaming- neutrino part 
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describes a phenomenological flow of neutrinos which freely stream out of the 
system^ 

2. A part of neutrinos produced by Eq. (|2-16p is assumed to leak out to be the 
streaming- neutrinos with a leakage rate Qj^ ak : 

V a (T»> s ) a p = Q l ^ k . (2-18) 

On the other hand, it is assumed that the remaining neutrinos constitute the 
trapped-neutrino part: 

Vp(T^t = Q a -Q 1 ^. (2-19) 

3. The trapped-neutrinos is combined with the fluid part as 

T aP = (T F ) Q/3 + (T^) af3 . (2-20) 
Then the equation for T a p is 

V Q T^ = -Qjf k . (2-21) 

We solve Eqs. (|2-18p and (|2-2ip . Note that the new equations only include 
the source term, Q^ ak , which is characterized by the leakage timescale ii ea k- 
Definition of Q£ 1S given in § 12.3.21 

The energy- momentum tensor of the fluid and trapped-neutrino parts (T a p) is 
treated as that of a perfect fluid, 

T a p = (j> + pe + P)u a up + Pg a p, (2-22) 

where p and u a are the rest mass density and the 4- velocity of the fluid. The 
specific internal energy (e) and the pressure (P) are the sum of the contributions 
from the baryons (free protons, free neutrons, a-particles, and heavy nuclei), leptons 
(electrons, positrons, and trapped-neutrinos) , and photons as, 

P = P B + Pe + P p h + P(u), (2-23) 

e = e B + e e + s ph + £(„) , (2-24) 

where subscripts , B\ 'e', 'ph', and '(z^)' denote the components of baryons, elec- 
trons and positrons, photons, and trapped-neutrinos (for v e , v e , and u x , see § 12.2.21) . 
respectively. Our treatment of the EOS is summarized in § 12.3.11 

The Euler equation (7™ V^T^ = — 7° (5„ ak ), and the energy equation (n Q V pT^ = 
—n a Q^ ak ) can be written explicitly, in terms of = T a ^n a n/s and ji = —T a P^i a np, 
as, 



= [~e h d t a + j k d^ k + ^S jk d ajk - , (2-25) 

dt{y/ie h ) + d k [^{e h v k + P(v k + (3 k ))~ 

= a ^7 (s^Kij - i k jid k In a + n^Gf*) , (2-26) 



*-* We note that Liebendorfer et alP^ developed a more sophisticate method in terms of the 
distribution functions of trapped and streaming neutrinos in the Newtonian framework. 
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where w = era* and v l = u' l /u t . 

The streaming-neutrino part, on the other hand, is written in the general form 

of 

{T u ' S ) a /3 = En a n p + F a np + F p n a + P af} , (2-27) 

where F a n a = P a pn a = 0. Then the evolution equations of streaming-neutrinos (E 
and Fj) are explicitly written as 

d t (^E) + d k [^(aF k - (3 k E)] = ^7 (<*P kl K kl - F k d k a - aQ l ^n a ) , (2-28) 

dt^Ff) + d k [^j{aP k - p k Fi\ 

= ^7 (~ma + F k dif3 k + ^P kl d ak i + aQf ak ) . (2-29) 

In order to close the system, we need an explicit expression of P a p (closure relation). 
In the current implementation, we adopt a simple form P a p = xEj a p with x = 1/3. 
We solve Eq. ()2T8p in a high resolution shock capturing schemed 

The closure relation employed in the current implementation is not very physical. 
Moreover we do not consider the so-called neutrino heating. To take into account the 
neutrino heating as well as the propagation of the streaming neutrinos accurately, a 
more sophisticated implementation of handling the neutrino transfer together with a 
better closure relation is required. However, such a study is beyond the scope of this 
paper. A more sophisticated implementation based on the moment methodj^'^ 
will be presented in the near future. 

2.2.2. Baryon and Lepton-number conservation equations 
The continuity equation for the baryon is 

V a (pu a ) = 0, (2-30) 

which can be written explicitly as 

dt(VlPw) + dki^pwv') = 0. (2-31) 
The conservation equations of the lepton fractions are written schematically as 

f = 7e, (2-32) 
dYr ' 7, e , (2-33) 



dt 

dt 
dY Ux 

dt 



7p„ (2-34) 
lu x , (2-35) 



where Y e , Y Ve , Y De , and Y Vx denote the fractions per baryon number for electrons, 
electron neutrinos (^ e ), electron anti-neutrinos ipe)-, and total of fi and r neutri- 
nos and anti-neutrinos (v x ), respectively. Note that only the trapped-neutrinos are 
responsible for these neutrino fractions. 
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Using the continuity equation for the baryon, we rewrite the conservation equa- 
tions of the lepton fractions in the following form as 

dt(VlpwY {L) ) + d k (^pwY {L) v k ) = V7«P7(L), (2-36) 

where Yas and 7^ are abbreviated expressions of the lepton fractions and the source 
terms. 

The source terms are given by 

(2-37) 
(2-38) 
(2-39) 
(2-40) 

where <-y local 's and 7 leak 's are rates of the local production and leakage for each species 
of neutrinos, respectively. As local reactions, we consider the electron capture, 
the positron capture, electron-positron pair annihilation, plasmon decay, and the 
Bremsstrahlung radiation of pair neutrinos (see § 12.3.21 for definitions and details). 
Because <-y local 's are characterized by the timescale of weak-interaction processes i wp , 
we follow the procedure proposed in Ref. 11) to stably solve the equations with a 
usual timestep (At ~ OA Ax) in an explicit manner (see Fig. [T]). 

1. At each timestep n, we first solve the conservation equation of the total lepton 
fraction (Y t = Y e + Y Ue - YpJ, 

^ = 7* = -(7^-7^), (2-41) 

instead of solving Eqs. lEEBZMZSD , together with Eqs. (|2l?ll . (12^25) - flZZgD , 
and (|2-35p . Note that the source term in Eq. (|2-4ip is characterized by the 
leakage timescale and can be solved explicitly. In this step, we assume that 
the /3-equilibrium condition is achieved, and the source terms are modified 
according to this assumption. After the time integration, the lepton fractions 
in the 'hypothetical' /3-equilibrium (Ye , Y„ e , and YpJ are calculated from the 
evolved Yj. 

2. We next solve the whole set of the equations (Eqs. {231]), fiERfy-^Rfy, and 
(|2-32p - (|2-35p ). In this step, we first calculate the maximum allowed values of the 
source terms, 7j,° C max an d 7p° C max' regarding that Y„ e and Yp e as the maximum 
allowed values of the neutrino fractions at the next timestep n + 1. Then the 
source terms are limited according to 



local 


= min 


local 

lu e J 


local 

7p e 


= min 


local 

!u e J 


^local 

Hue 


= min 


^-jlocal 
Hv e 


^-jlocal 
Vp e 


= min 
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local 

/V e ,max 

local 

7p e ,max 



(2-42) 
(2-43) 
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Fig. 1. A schematic picture of the leakage scheme. See text for details. 



These limiter procedures enable us to solve the equations in an explicit manner. 
3. After the evolution, following conditions are checked, 



(2-46) 
(2-47) 



where fi p , fj, n , fj, e , fj, Ue , and /ip e are the chemical potentials of protons, neutrons, 
electrons, electron neutrinos, and electron anti-neutrinos, respectively. If both 
conditions are satisfied, the values of the lepton fractions at the timestep n + 1 
are reset to be those in the /3-equilibrium value; Y<P , Yj? e , and Yp e . 

2.3. Microphysics 
2.3.1. Equation of state 

We employ two versions of Shen's EOS. One is a purely nucleonic EOS,® which 
is adopted both in stellar core collapse simulations (see § [3]) and in BNS merger 
simulations (see §0|). In the BNS merger simulation, we also adopt an EOS in which 
effects of A hyperons are taken into accounPS (hereafter, referred to as Hyp-EOS). 
These EOS are tabulated in terms of the rest-mass density (p), temperature (T), 
and Y e or Y\. 

Shen-EOS, derived from a relativistic mean-field theory,^ 1 is a stiff one among 
many other EOS, giving a large maximum gravitational mass of zero-temperature 
spherical neutron stars M max ~ 2.2M Q . By contrast, M max ~ 1.8Mq for Hyp-EOS 
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because the appearance of hyperons softens the EOS. The latest discovery of a high- 
mass neutron star with mass 1.97 ± 0.04M<J^ suggests that stiff EOS are favored, 
and Shen-EOS satisfies this requirement whereas Hyp-EOS does not. However, we 
consider that Hyp-EOS is a viable candidate for the neutron-star EOS except for a 
very high density that an only high- mass neutron star of M < M max has. We note 
that the neutron star in the BNS and a hypermassive neutron star (HMNS) studied 
here do not have the extremely high density (except for the HMNS just before the 
collapse to a BH). 

The thermodynamical quantities of a dense matter at various sets of (p, Y p , T) 
are calculated to construct the numerical data table for simulations. Here Y p is the 
total proton fraction per baryon number. The original table covers a range of density 
10 51 -10 15,4 g/cm 3 , proton fraction 0.0-0.56, and temperature 0-100 MeV, which are 
sufficient parameter ranges for supernova simulations. The original table has been 
extended to a higher density (10 51 -10 17 g/cm 3 jSSi^B anc i a higher temperature 
(0-400 MeV)^>^ for following the BH formation. 

It should be noted that the causality is guaranteed to be satisfied in this frame- 
work, whereas the sound velocity sometimes exceeds the speed of the light in the 
non-relativistic framework, e.g., in the EOS by Lattimer and SwestyPS This is one 
of the benefits of the relativistic EOS. 

To consistently calculate the pressure and the internal energy of electrons and 
positrons, the charge neutrality condition Y p = Y e has to be solved to determine the 
electron chemical potential p e for a given set of p and T in the EOS table. Namely, 
it is required to solve the equation 

n e (p e ,T) = n_ - n + = — (2-48) 

in terms of p e for given values of p, T, and Y e (= Y p ). Here, m u = 931.49432 MeV is 
the atomic mass unit, and n_ and n+ are the total number densities (i.e., including 
electron-positron pairs) of electrons and positrons, respectively. Then, assuming that 
electrons and positrons obey the Fermi-Dirac distribution, the number density, the 
pressure, and the internal energy density of electrons and positrons are calculated.^ 
The pressure and the specific internal energy density of photons are given by 

P ph = ^, e ph = ^ll, (2-49) 
d p 

where a r = (ir 2 kg) / (15c 3 h 3 ) is the radiation constant. 

In our leakage scheme, the trapped-neutrinos are assumed to interact sufficiently 
frequently with the matter that be thermalized, and hence, they are described as ideal 
Fermi gases with the matter temperature. From the numerically evolved trapped- 
neutrino fractions Yr v \ (u e , u e , and v x are abbreviated by (u)), the chemical potentials 
of the trapped-neutrinos (fJ>r u <\) are calculated by solving 

Y( V ) = ^yn (u) (p {u) ,T), (2-50) 

where nt v \ is the number density of the trapped-neutrinos. Then the pressure and 
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the internal energy of the trapped-neutrinos are calculated in the same manner as 
for electrons, using p^ and the matter temperature. 

In the high-resolution shock-capturing scheme for hydrodynamics, we in general 
need to evaluate the sound velocity c s , 



(2-51) 
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Here, the derivatives of the pressure are calculated by 



dP 
dp 

dP 
d7 



E 



E 



dp 
dT 



r 



dPi 
dT 




dej_ 
dp 



-i 




(2-52) 
(2-53) 



where the sum is taken over B,e,ph and (y). 
2.3.2. Weak-interaction and leakage rates 

The leakage rates are phenomenologically defined by^ 



'^(if) - .local 



lleak 

■J) u a 



i{y) 



_ leak 

7 M 



(l_ e -^) )7 difF + e - 



7 M 



(2-54) 
(2-55) 



where 77^ is the optical depth of neutrinos and 6 is a parameter which is typically set 
to be b~ 1 = 2/sEjj Note that <3'°^ k should be regarded as the emissivity of neutrinos 
measured in the fluid rest frame so that we set Qa ak = Q l ^u a ^'^ 
The optical depth is calculated by2D'E3 



T(u) = mln 

r (iy ) = min 



_x y z 



(2-56) 
(2-57) 



for axisymmetric (see § [3]) and three-dimensional (see § UJ simulations, respectively. 
Here r^, r^, r^, and are the optical depths along x, y, z, and a 'quasi-radial' 



directions from each grid point, respectively. We calculate, for example, t? ^ by**) 

T ( ^(tu,z) = J B M (ro,2:) 2 f z (ro,2), (2-58) 

ut 

k(w,z')oIz', (2-59) 



f 2 (tu, z) 



*' We again used an abbreviation {v) for u e , u e , and v x 
* T (v) ^ r (V) are calculated in a similar manner. 
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where z out denotes the outer boundary in the z-direction. R (= K(„)/E?s) is an 
'opacity' in which the neutrino-energy dependence is factored out (see Appendix C). 
The neutrino energy is determined by 

E {v) = (1 - e - T M /c )£g + e - T ^ /c E$f, (2-60) 

where we set the parameter as c = 5, which implies that it takes about three collisions 
to thermalize a neutrinosP^ Note that T(„) depends on the neutrino energy Ef v \ 
and we solve this equation by the Newton- Raphson method. Effi and E^ al are 
the average (thermalized) diffusion and local-production energy, which are given 
respectively by, 

diff _ F 3 (p M /k B T) 
E, u) - k B l — n > (2-blJ 

/olocal 

P 7 ( „) 

where Fk{x) is the Fermi-Dirac integral. 

As the local production reactions of neutrinos, we consider the electron and 
positron captures (7^ and 7p ( f),' the electron-positron pair annihilation (7^^ for 
electron- type neutrinos and j^d x f° r other types) JSSJ the plasmon decays (t£J^ and 
7^'p^)|2D and the Bremsstrahlung processes (7^ e r p™ s and 7^Jp™ s )-^ Then, the local 
reaction rates for the neutrino fractions are 

local ec , pair plas Brans (2 .qo\ 
local pc pair pla_s Br_ems ^.g^ 
local = 4 / pair plas BremsN ^.g^ 

Similarly, the local neutrino energy emission rate Q\° u c f is given by 

Ql° e cal = Qll + (Qlie + Qtl + Qu T JD . ( 2 ' 66 ) 
Ql° cal = Q£ + {Ql% + Ql l 2 + Q* r JT s ) > (2-67) 

Ql°; al = 4 (Q^t + Qg£ + Q^ ms ) • (2-68) 

The explicit forms of the local rates in Eqs. (|2-63p - (|2-68p are summarized in Appen- 
dices A and B (see also, Ref. fTTp ). 

We follow the recent work by Rosswog and Liebendorfer 22 ' for the diffusive 
neutrino emission rates 7^ and Qgf in Eqs. (|2 T MD and (12^551) . The explicit forms 

of 7^ and Q^f are described in Appendix C (see also, Ref. [TTj) ). 

2.4. Recover of (p, Y e /Y h T) 

The quantities numerically evolved in the hydrodynamic equations are the con- 
served quantities: ^ypw, ^ypwYL, yFfji, and 1 /7"e/ l . The argument variables, (p, 
(Y e or Yi), T), of the EOS table, together with w = otu 1 = yl + j^U{Uj, should be 
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calculated from the conserved quantities at each timestep. Note that ^7 is readily 
given by numerical evolution of Einstein's equations. Also, the lepton fractions (Yl) 
can be calculated directly from the conserved quantities. 

2.4.1. Non-/3-equilibrium case 

In the case that the /3-equilibrium condition is not satisfied, the argument quan- 
tities (p, Y e , T) can be reconstructed from the conserved quantities in the following 
straightforward manner. 

1. Give a trial value of w, referred to as w. Then, one obtains a trial value of the 
rest mass density p. 

2. A trial value of the temperature, T, can be obtained from the numerically 
evolved value of e^, by solving the following equation: 

e h -Ys e h,(.)(P,Y M f) = e htEOS (p,Y e ,f). (2-69) 
M 

Here, eh,EOs{p~,Y e ,T) = pw 2 hEos{f>,Y e ,T) — P{p,Y e ,T) should be evaluated 
from the EOS table which does not include the contributions of trapped-neutrinos. 
In the left-hand-side, ZhSv) is the trapped-neutrino part. Note that one dimen- 
sional search over the EOS table is required to obtain T. 

3. The next trial value of w is given by 

w = ,/l + 7 « (^—) i^J—). (2-70) 
V \pwh E osJ \pwh E osJ 

4. Repeat the procedures (l)-(3) until a required degree of convergence is achieved. 
Convergent solutions of the temperature and w are obtained typically within 
10 iterations. 

2.4.2. The /3-equilibrium case 

In the case that the /3-equilibrium condition is satisfied, on the other hand, we 
may reconstruct the argument quantities (p, Y e , T) from the conserved quantities and 
Yi, under the assumption of the /3-equilibrium. In this case, two-dimensional recover 
(Yi,eh) ==> (Y e ,T) would be required for a given value of w. In this case, there 
may be more than one combination of (Y e , T) which gives the same values of Y\ and 
eh- Therefore, we have to adopt a different method to recover (p, Y e ,T). Under the 
assumption of the /3-equilibrium, the electron fraction is related to the total lepton 
fraction: Y e = Y e (p,Yi,T). Using this relation, the EOS table can be rewritten in 
terms of the argument variables of (p, Yi, T). Then, the similar strategy as in the 
non-/3-equilibrium case can be adopted. Namely, 

1. Give a trial value w. Then one obtains a trial value of the rest mass density. 

2. A trial value of the temperature can be obtained by solving 

e h = 4 os (p,Y l ,f) (2-71) 

with one dimensional search over the EOS table. Here e^ os should be evalu- 
ated from the /3-equilibrium EOS table, which contains the trapped-neutrino 
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contributions. 

3. The next trial value of w is given in the same way. 

4. Repeat the procedures (l)-(3) until a required degree of convergence is achieved. 
The electron fraction is given as Y e = Y e (p,Yi,T) in the /3-equilibrium EOS 
table. 

In the case of a simplified or analytic EOS, the Newton-Raphson method may 
be applied to recover the primitive variables. In the case of a tabulated EOS, by con- 
trast, the Newton-Raphson method may not be a good approach because it requires 
derivatives of thermodynamical quantities which in general cannot be calculated 
precisely from a tabulated EOS by the finite differentiating method. 

§3. Gravitational collapse of massive stellar core 

The observational associations (for a review, see Ref. [60]) ) between LGRBs and 
supernovae has provided the strong support to a scenario, so-called collapsar model, 
in which LGRBs are assumed to be driven in the collapse of a massive stellar core to 
a BrljSDJSZJ in the collapsar model, a central core of a massive star is required to be 
rotating rapidly enough that a massive accretion disk can be formed around a BH. 

Because the observed supernovae associated with LGRBs are Type Ib/c and the 
relativistic jets have to reach the stellar surface,^ the progenitors should have lost 
their hydrogen (and helium) envelopes before the onset of the stellar core collapse; 
otherwise a peculiar evolution path is required. Due to these reasons, the progenitors 
of LGRBs are now believed to be rotating massive Wolf-Rayet (WR) stars. However, 
ordinary WR stars are known to be accompanied by strong stellar winds driven by 
the radiation pressure which cause a rapid spin-down of the stellar core. Here, a 
serious problem concerning the collapsar model is that according to stellar evolution 
calculations, it is very difficult to produce pre-collapse cores which satisfy both the 
requirement of the collapsar model and the association of Type Ib/c supernova, if 
magnetic torques and standard mass-loss rates are taken into account.^ 1 

To resolve the above dilemma, several models have been proposed (see Ref. |65|) 
for a review). All of the proposed progenitor models of LGRBs are anomalous in 
the sense that they are different from the progenitors of ordinary supernovae (see 
Ref. [T2]) for a discussion). Qualitatively speaking, LGRB progenitor cores may be 
modeled by a rapidly rotating, higher-entropy core, regardless of their formation 
processes. Based on this assumption, we performed simulations of a massive stellar 
core with higher values of entropy collapsing to a BH. In this section, we report our 
latest results of fully general relativistic simulations for the collapse of a rotating, 
higher-entropy cores, performed taking into account detailed microphysics. 

3.1. Initial models and grid setting 

As a representative model of a high entropy core, we adopt a presupernova core 
of IOOMq model calculated by Umeda and Nomotd^ 1 (hereafter denoted by UN100). 
The model has an iron core of a large mass M COTC ~ 3.2M & and radius -R CO rc ~ 2500 
km with the central density and temperature of p c ~ 10 9 ' 5 g/cm 3 and T c ~ 10 10 K. 
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<P C < 0.02 


< <P C < 0.044 


< $c < 0.09 


<<P C < 0.2 


<P C > 0.2 


Axo (km) 


4.0 


2.0 


1.0 


0.5 


0.25 


S 


0.0075 


0.007 


0.0065 


0.006 


0.0055 


N 


332 


428 


542 


668 


812 


L (km) 


5840 


5370 


5000 


4450 


3860 


Axo (km) 


5.5 


2.85 


1.45 


0.7 


0.35 


8 


0.0075 


0.007 


0.0065 


0.006 


0.0055 


N 


294 


380 


484 


610 


752 


L (km) 


5860 


5360 


4980 


4370 


3870 



Table I. Summary of the regridding procedure. The values of the minimum grid spacing Axq (in 
units of km), the non-uniform-grid factor 8, and the grid number N for each range of $ c — 1 — a c 
are listed for the finer and the coarser (lower table) resolutions. 



The central value of entropy per baryon is s ~ which is much larger than that 
of an ordinary presupernova core for which s < Iks- 

Because the model UN100 is non-rotating, we add rotational profiles according 

R 2 

Q(w) = i? 2 — ^Jcut, (3-1) 
Rq + zu z 



where w = \J x 2 + y 2 , Oq, and Ro are parameters which control the magnitude and 
degree of differential rotation. The cut-off factor is introduced by a practical 
reason for the numerical simulation: If the specific angular momentum in the outer 
region of the core is too large, the matter escapes from the computational domain. 
To avoid this, the rotational velocity has to be suppressed in the outer region. 

We fix the central angular velocity as Qq = 1.2 rad/s and consider two values of 
-Ro; a rigid rotation model {Rq = oo, referred to as UNlOO-rigid) and a differential 
rotation (Rq = R COTe , referred to as UNlOO-diff) model. We note that the imposed 
rotation is moderately large (not rapid) because Qq is much smaller than the Kepler 
value (M core /R^ ore y/ 2 rj 5.2 rad/s. 

We assume axial and equatorial symmetries of the spacetime and the so-called 
Cartoon method^'^ is adopted for integrating Einstein's equations. In the cur- 
rent implementation, we use a fourth order Lagrange interpolation scheme, which is 
necessary in the Cartoon method. 

In numerical simulations, we adopt a nonuniform grid, in which the grid spacing 
is increased according to the rule 

dxj + \ = (1 + 5)dxj, alzi + i = (1 + 5)dzi, (3 - 2) 

where dxj = Xj+i — Xj, dz\ = z\+\ — zi, and 5 is a constant. In addition, a regridding 
technique^ >EQ j s adopted to assign a sufficiently large number of grid points inside 
the collapsing core, saving the CPU time efficiently. The regridding is carried out 
whenever the characteristic radius of the collapsing core, defined byZ® <P C = 1 - 
a c > 0) where a c is the central value of the lapse function, decreases by a factor 
of ~ 2, and we set an infalling boundary condition at the outer boundary. 

All the quantities on the new grid are calculated using a fifth-order Lagrange 
interpolation. However, for the fluid quantities such as p and h, the fifth-order 
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interpolation could fail because the interpolation may give negative values of p and 
h — 1. In case we have p < or h < 1, we adopt the linear interpolation to calculate 
the quantities on the new grid, based on the prescription proposed by Ref. I42|) . In 
each regridding, we solve the Hamiltonian constraint equation numerically. 

To check the convergence of numerical results, simulations are performed in two 
different grid resolutions. Table |T] summarizes the regridding parameters (N and 
L are the number of the grid points and the computational domain, respectively) 
of each level of the regridding procedure for finer (upper) and coarser (lower) res- 
olutions. The numerical results in both grid resolutions agree well except for the 
formation time of a BH and the stochastic behavior due to connective and turbulent 
motions. 

3.2. Dynamical features 

Figure 12(a) shows the evolution path (red curve) of central values of the rest- 
mass density and the temperature in the p-T plane for UNlOO-rigid. As in the core 
collapse of an ordinary supernova for which the central value of entropy per baryon is 
s/ks ~ 1, gravitational collapse is triggered by the electron capture and the photo- 
dissociation of heavy nuclei. Because of the higher value of the entropy per baryon 
(s/ks ~ 4), the photo-dissociation is mainly responsible to the destabilization. Note 
that a substantial amount of heavy nuclei are resolved into heliums by the photo- 
dissociation (see Fig. [2]): The fraction of heavy nuclei in mass is ~ 0.4 and 0.2 for 
p c = 10 11 and 10 12 g/cm 3 , respectively. Then the collapse in the early phase proceeds 
in a homologous manner. As the collapse proceeds temperature increases, heliums 
are resolved into free nucleons (p, n). 

Figure [2(b) shows the evolution path (red curve) of central values of the rest- 
mass density and Y e in the p-Y e plane for UNlOO-rigid. Because the temperature 
for the present models is higher than that for the ordinary supernova, the electron 
capture on the free proton is enhanced due to the larger value of the free proton 
fraction, and hence, the electron fraction for UN100 is by ~ 0.1 smaller than that 
for the ordinary supernova in the collapse phase (compare the red and green curves 
in Fig. [2(b)). 

The time evolution of the central values of the rest-mass density, electron frac- 
tion, temperature, and the lapse function for models UNlOO-rigid and UNlOO-diff is 
shown in Figure [31 As in the collapse of the ordinary supernova core, the collapsing 
core experiences a bounce when the central density reaches the nuclear density p nuc 
above which the pressure increases drastically due to the repulsive nuclear force, 
and then, shock waves are formed and launched. Because the electron fraction at 
the bounce is small as F e ~ 0.17 and hence the core is neutron rich, the nuclear 
force starts playing a role at relatively low density, p ~ 10 14 g/cm 3 , in Shen-EOS. 
After the bounce, a HMNS, which is supported by a significant rotation and thermal 
pressure, is formed. 

The shock wave formed at the core bounce propagates outward but eventually 
stalls at r ~ 100 km due the neutrino cooling and photodissociaion of heavy nuclei 
contained in the infalling matter (see the top panels of Figs. U] and [5]). Then, a 
standing accretion shock is formed. As in the case of the ordinary supernova, con- 
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(b) 




Density [g/cm ] 



Fig. 2. (a) Evolution paths of the central values of the rest-mass density and the temperature in 
the p-T plane. The red curve shows the evolution path of UNlOO-rigid. The black solid curve 
shows the boundary at which the condition P e — P gas is satisfied (P e > P gas for the higher 
density side). The two blue dashed curves denote the values of (p,T) with which 56 Fe or 4 He 
will be half by mass due to the photo-dissociation. An evolution path for an ordinary supernova 
core^ is shown together for comparison (solid green curve), (b) Evolution paths of the central 
values of the rest-mass density and the electron fraction in the p-Ye plane. 
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Time [ms] 

Fig. 3. Time evolution of the central values of rest-mass density, electron fraction, temperature, 
and the lapse function for UNlOO-rigid (blue curves) and UNlOO-diff (red curves). The results 
for the coarser grid resolution are shown together (thin curves) . The collapsing core experiences 
a core bounce at t ~ 510 ms and collapse to a BH at t ~ 1810 ms. 

vection is activated between the HMNS and the standing shock. However, it is not 
strong enough to push the standing shock outward. The convection is stronger for 
the differentially rotating model, which rotates more slowly than the rigidly rotating 
model. This is likely to be due to the stabilizing effect of the epicyclic modes which 
is stronger in UNlOO-rigid. 121 As the matter accretion proceeds, the central density 
and temperature increase gradually, and eventually, the HMNS collapses to a BH. 
The formation time of the BH depends on the grid resolution (compare the thin and 
thick curves in Fig. [3|) in particular for UNlOO-rigid. This is because the HMNS is 
close to the marginally stable configuration, and hence, a small thermodynamical 
change results in a significant change in p. In general, for a finer grid resolution, 
the lifetime of HMNS increases. The longer lifetime for higher-resolution runs is a 
often-seen feature, because the numerical dissipation is less severe for high resolution. 
Note that the rotational profiles in the central region of UNlOO-rigid and UNlOO-diff 
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X [km] 

Fig. 4. Contour profiles of the rest-mass density in the x-z plane at t — 570 ms, (top left), 645 ms 
(top right), 890 ms (bottom left), and 1150 ms (bottom right) for UNlOO-difT. 

are very similar, and their evolution process agrees well with each other soon after 
the core bounce. However, their evolution paths deviate as the matter in the outer 
region falls. 

Dynamics of the system for the models UNlOO-diff and UNlOO-rigid in the later 
phase of the accretion onto the HMNS is qualitatively different. Figure H] shows 
contour plots of the rest-mass density in the x-z plane at selected time slices for 
UNlOO-diff until the HMNS collapses to a BH. Because of the accretion of the matter, 
the shock front of the standing accretion shock gradually recedes. Note that the 
shape of the shock wave is deformed by the rotation to be spheroidal. As we shall 
see below, this shows a remarkable contrast with the case of UNlOO-rigid where the 
shock wave is deformed to be a torus-like shape. When the shock wave stalls, negative 
gradients of the entropy per baryon and of the total-lepton (electron) fraction appear 
because neutrinos carry away both the energy and the lepton number, as in the 
collapse of an ordinary presupernova core (see the bottom left panel of Fig. The 
HMNS finally collapses to a BH due to the mass accretion, and a geometrically thin 
disk is formed around the BH. The system shows no violent time variability. 
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Fig. 5. Contour profiles of the rest-mass density in the x-z plane at t — 617 ms, (top left), 863 ms 
(top right), 1256 ms (middle left), 1437 ms (middle right), 1822 ms (bottom left), and 2225 ms 
(bottom right) for UNlOO-rigid. In the middle left panel, outflow velocity vectors larger than 
0.15c are plotted together (red arrows). Note that the scale of the figure is different from that 
of Fig. [J] In the bottom two panels, a BH is formed at the center. 
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Fig. 6. Contour profiles of the electron fraction in the x-z plane at the same timeslices as Fig. [S] 

Figures El El and El respectively, show contour plots of the rest-mass density, 
the electron fraction, the entropy per baryon, and the temperature in the x-z plane 
at selected time slices for UNlOO-rigid. As in the model UNlOO-diff, the shock wave 
formed after the core bounce stalls at r ~ 100 km (the top left panel of Figs. [5]- 
[8]). Due to the faster rotation of the outer region than in UNlOO-diff, the shock 
wave is deformed to be a torus- like configuration (the top right panel of Figs. [SHED- 
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X [km] 

Fig. 7. Contour profiles of the entropy per baryon in the x-z plane at the same timeslices as Fig. 



The formation of this torus-shaped shock is the key ingredient which characterizes 
the dynamics of UNlOO-rigid. At the shock, the kinetic energy associated with the 
motion perpendicular to the shock surface is dissipated but that associated with 
the parallel component is preserved. In the model UNlOO-rigid, the shock front is 
highly deformed, and thus, the amount of the kinetic energy dissipated at the shock 
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Fig. 8. Contour profiles of the temperature in the x-z plane at the same timeslices as Fig. [S] but 
with a different (zooming) scale. 

is not as large as that in UNlOO-diff. This implies that the infalling materials are 
eventually accumulated in the central region and their kinetic energy is dissipated 
at the surface of the HMNS. Figure O which displays the velocity field in the x-z 
plane at a time slice, clearly shows this mechanism. During this process, oscillations 
of the HMNS are excited as the infalling matter hits it. Also, the shock waves gain 
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X [km] 

Fig. 9. Velocity vector fields outside the shock surface (blue arrows) and inside the shock (red 
arrows) together with a contour profile of the rest-mass density in the x-z plane at t = 645 ms. 
Only the velocity component perpendicular to the shock surface is dissipated at the shock. As 
a result, the infalling matter is accumulated in the central HMNS. 

the thermal energy via PdV work and propagate outward. 

Due to the accumulation of the matter onto the HMNS and the resulting shock 
heating, the thermal energy is stored in the polar region of the HMNS, increasing 
the gas pressure, P gas , there. On the other hand, the ram pressure, P r am> of the 
infalling matter decreases with the elapse of the time because its density decreases. 
When the condition, P ram < P gas , is realized, outflows are launched from the polar 
surface of the HMNS, forming shocks (see the middle left panel of Figs. [5H8|). It can 
be seen that the entropy around the rotational axis is significantly enhanced due to 
the shock heating associated with the outflows (e.g., see the middle right panel of 

Fig. ED. 

The outflows eventually lose the driving power by the neutrino cooling and mat- 
ter again turns to fall onto the polar region of the HMNS. Due to the continuous 
mass accretion, the HMNS eventually collapses to a BH surrounded by a geometri- 
cally thick torus (see the bottom left panel of Figs. [5HH]). Note that in the present 
leakage scheme, neutrino heating is not taken into account and exploring the fate 
of the thermally driven outflows in the presence of the neutrino heating is an in- 
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teresting subject. We plan to pursue this issue using a code based on the moment 
formalism®}'^ in the near future. 

We found, as another novel feature of dynamics, that the BH-torus system shows 
a time variability (see the bottom right panel of Figs. [SHSD ■ This is reflected in the 
neutrino luminosities as we shall show in § 13.31 Such a time variability has not been 
seen in UNlOO-diff. Reason for this will be explained as follows. 

First, the infall timescale of the matter in the torus into the BH is longer for 
UNlOO-rigid due to the rapider rotation (in the outer region). Also, the neutrino 
cooling timescale for the torus will be longer for UNlOO-rigid because of the larger 
optical depth due to the higher density and temperature. Furthermore, the heating 
rate due to the mass accretion in the central region is larger for UNlOO-rigid due 
to the mass accumulation mechanism. Due to these reasons, the energy deposition 
by the accretion cannot be valanced by these cooling mechanisms: Q& cc > Q^^u + 
Q~. Then the torus will expand lowering the optical depth which results in the 
enhancement of the neutrino cooling rate Q~ . Because of the strong dependence 
of neutrino opacities and cooling rate on the temperature, a slight change in the 
shock configuration may result in a huge loss of the thermal energy by neutrino 
emission. Here, note that the shock heated matter is partially supported by the 
pressure gradient due to the moderate (not very rapid) rotation. Therefore, if some 
materials lose their thermal energy, they will drop into the BH like an avalanche. The 
modulation of the shock configuration which triggers the above dropping appears to 
come from the Kelvin-Helmholtz instability, developed at the interface between the 
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torus and the accumulating flows. 

In the model UNlOO-diff, by contrast, Q^ cc is smaller due to the absence of the 
accumulation mechanism, and QiLyj and Q~ are larger due to the slower rotation. 
As a result, the above energy balance will be satisfied without expansion of the disk, 
and thus, there is no violent time variability. 

It is remarkable that the above qualitative differences in dynamics between 
UNlOO-diff and UNlOO-rigid stem from a small difference in the initial angular veloc- 
ity profile in the outer region. Figure [10] compares profiles of the rotational angular 
velocity along the equator just before the BH formation for UMlOO-diff and UM100- 
rigid. The rotational profiles of the HMNS (r < 40 km) are similar in the central 
region, and in the outer region, the difference at most by a factor of 2-3. This result 
shows that the final outcome depends strongly on the rotational profile of progenitor 
stars. 

3.3. Neutrino Luminosity and Gravitational Waves 

Figures [TlTa) and flHb) plot the time evolution of neutrino luminosities for 
UNlOO-diff. In the prebounce phase, electron neutrinos are dominantly emitted 
and the emissivity of electron anti-neutrinos is much smaller. This is because the 
electrons are (mildly) degenerate blocking the inverse /3-decay and also the positron 
fraction, which is responsible for the anti-neutrino emission, is small. Soon after the 
core bounce, the so-called neutrino burst occurs at the time that the shock wave 
passes through the neutrino-sphere, as in the collapse of an ordinary supernova core. 

After the neutrino burst, the emission of electron anti-neutrinos is enhanced 
and their luminosity becomes larger that of electron neutrinos. This property is 
different from that in the ordinary supernova,^® and explained as follows. During 
the post neutrino burst phase, a large number of positrons are produced because the 
degeneracy parameter becomes low as r/ e ~ 1 due to the high temperature of T > 20 
MeV, which is higher than the temperature in the ordinary supernova T ~ 5 MeV.^ 
Then, because the neutron fraction X n is much larger than the proton fraction X p , 
the positron capture on neutrons occurs more efficiently, and hence, the electron 
anti-neutrino luminosity becomes larger than the electron neutrino luminosity. This 
dominant emission of electron anti-neutrinos are also found for a HMNS formed after 
the BNS merger (see § 0]) and in a BH-torus system formed in the collapse of a more 
massive core^ with s = 8/cb. Both systems have a higher temperature and a lower 
electron fraction than the collapse of the ordinary supernova core, as in the present 
case. 

The luminosity of \x and r neutrinos is smaller than that of electron neutrinos and 
anti-neutrinos. This is simply due to the absence of the neutrino production channel 
mediated by the charged weak current. At later phases in the HMNS evolution (800 
ms < t < 1100 ms), the electron neutrino and anti- neutrino luminosities show weak 
time variability. This is due to the convective activities that occur near the neutrino 
sphere. The fi and r neutrino luminosity does not show the variability because they 
are mainly emitted by the hot central regions that do not suffer from the convection. 

Soon after the BH formation at t ~ 1090 ms, neutrino luminosities decrease 
drastically because the main neutrino-emission region is swallowed into the BH. After 
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Fig. 11. Time evolution of neutrino luminosities for UNlOO-diff (a) before the BH formation and 
(b) after the BH formation. The red, green, and blue thick curves correspond to the luminosities 
of u e , i> e , total of fj, and r pair neutrinos, respectively. The thin blue curve in the upper panel 
shows the luminosity of individual h/t neutrino (namely the quarter of thick blue curve). The 
thick black curve in the lower panel shows the total neutrino luminosity. 
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Fig. 12. Time evolution of neutrino luminosities for UNlOO-rigid (a) before the BH formation and 
(b) after the BH formation. Meanings of all curves are the same as those of Fig. [11] 
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that, the geometrically thin accretion disk emits ~ 10 51 -10 52 ergs/s by neutrinos in 
its early evolution phase with the duration ~ 100 ms and the luminosities decrease 
monotonically in time. We do not find any enhancement of the neutrino luminosities 
after the BH formation in our simulation time. In this phase, electron neutrinos 
are dominantly emitted because the disk is at a lower temperature of T < a few 
MeV, and hence, there are less positrons. Also, the number of the target neutrons 
are smaller because the disk is composed mainly of proton-rich matter in the outer 
neutrino emission region. 

Figures fT2Ta) and [T2T b) show the time evolution of neutrino luminosities for 
UNlOO-rigid. The features of the evolution is similar to those for UNlOO-diff be- 
fore the neutrino burst (t < 700 ms). After that time, the luminosities of electron 
neutrinos and anti-neutrinos gradually decrease. This is because the optical depth 
(diffusion time) gets larger (longer) as the torus grows. The luminosities of electron 
neutrinos and anti-neutrinos show only weak time variability, reflecting the weaker 
convective activity in UNlOO-rigid. 

At t « 800 ms, the total luminosity of \x and r neutrinos (L v + Lp M + L Vr + L Dt ) 
becomes larger than the electron neutrino and anti-neutrino luminosities. This is 
partly due to a very high temperature of the neutrino sphere which enhances pair 
neutrino production processes, as well as a smaller optical depth along the rotational 
axis: Thermal neutrinos from the hot HMNS will be almost directly seen due to the 
low density along the rotational axis. Indeed, after t > 1300 ms, all flavor of neutrinos 
and anti-neutrinos are almost equally emitted, indicating the dominant emission of 
thermal neutrinos from the very hot HMNS in the neutrino luminosity. This feature 
is not seen in the BNS merger (see § and the collapse of the more massive stellar 
core JI2J ^his j s because the continuous mass accretion is absent in the BNS merger 
and the HMNS is quickly collapses to a BH in the collapse of the more massive 
stellar core. At the final phase in the fallback collapse of an ordinary core (a failed 
supernova),^ 1 an enhancement of the emission of v x is also seen. 

Note that the above result implies that observational signals of neutrinos could 
depend on the viewing angle: If we would see the system from the direction along the 
rotational axis, we might see a brighter emission of neutrinos with a higher average 
energy from the hot HMNS, while we would see neutrino emissions from the torus 
if we see the system from the direction along the equator. In the leakage scheme 
adopted in this paper, unfortunately, we cannot investigate such an angle dependence 
of neutrino luminosities. 

The neutrino luminosities show a precipitation when the BH is formed at t ~ 
1805 ms, as in the case of UNlOO-diff. The total neutrino luminosity emitted from 
the torus around the BH amounts to L„,tot ~ 10 51 ergs/s. A remarkable property in 
UNlOO-rigid is that this luminosity is maintained for > 1 s. In addition, by contrast 
with the case of UNlOO-diff, the neutrino luminosities show a violent time variability. 
Such a long-term high luminosity and a time variability may be associated with the 
time variability that LGRBs show. 
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Fig. 13. The spectrum of the characteristic gravitational-wave strain for UNlOO-rigid. The noise 
amplitudes of Advanced LIGO for a version in which no signal recycling mirror is used (NO 
SRM) and a design with a narrow-band tuning at 1kHz (High Freq) are shown together. 



Figure [T3] plots the spectra of the characteristic gravitational- wave strain,^ 

(3-3) 
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is the energy power spectra of the gravitational radiation. Ai (/) is the Fourier 
transform of A2, 

Mf) = J A 2 (t)e 2 ^ t dt (3-5) 
with A2 being the +- mode of gravitational waves with I = 2 and m = 0, 
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where Ijj denotes a quadrupole moment,^ 1 I{j its second time derivative, and t ret a 
retarded time. Because the strain of gravitational waves has a broad-band spectrum, 
they appear to be emitted primarily by a long-term stochastic motion of the infalling 
material and of the matter in the HMNS, as in the collapse of ordinary coresP^ 1 
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Note that /i c har includes only gravitational waves from the matter contribution and 
not from the anisotropic neutrino emissions. We also show the noise amplitudes of 
Advanced LIGO for a version in which no signal recycling mirror is used (NO SRM) 
and of a version of a specially designed narrow-band Advanced LIGO (High Freq) 
together P3 

The effective amplitude of gravitational waves observed in the most optimistic 
direction is /i c har ~ 10 -20 for an hypothetical event at a distance of 10 kpc, which is 
larger than that for the ordinary supernovaP^ The frequency at the peak amplitude 
is 1-2 kHz. The reason for this larger gravitational- wave amplitude is that the 
HMNS in the present model has a larger mass and the gravitational-wave signal is 
accumulated during the long-term convective activity, for which the duration > 1 s 
is longer than that for the ordinary supernova. Figure [TBI shows that with NO-SRM 
Advanced LIGO the signal of gravitational waves will be detected with a signal-to- 
noise ratio S/N ~ 10 for D = 10 kpc. For High- Freq Advanced LIGO, the detection 
may be done with S/N ~ 100 for D = 10 kpc. Even for D = 50 kpc (the distance to 
the Large Magellanic Cloud), the detection will be possible with S/N > 20. 

For an event of D < 50 kpc, a large number of neutrinos will be also detected by 
water-Cherenkov neutrino detectors such as Super-Kamiokande and future Hyper- 
Kamiokande (see also a discussion in § I4.2|) . 

3.4. Possible association to Gamma-ray bursts 

Before closing this section, we give an order estimate of the energy deposition 
rates (E u &) by the neutrino pair-annihilation which is one of the possible processes 
to drive relativistic jets required to produce LGRBs. Note that the energy from 
the neutrino pair annihilation should be deposited in a baryon-poor region in order 
to generate highly relativistic outflows. The funnel region near the rotational axis 
above the torus formed in UNlOO-rigid is a promising place for this purpose. 

According to the estimate in Ref.l76p. the deposition rate in the BH-torus system 

would be proportional to M 9 / 4 M B y //2 . In this estimation, the neutrino luminosity 

is assumed to be originated from a viscous heating. In our present simulations, 

the neutrino luminosity is determined by the infalling rate of the material which 

experiences the shock heating at the surface of the torus to increase the thermal 

energy of the torus. However, the dependence of the pair-annihilation rate on the 

mass infall rate M is essentially the same for thick torus phase. Figure [TH shows the 

irreducible mass of the BH (the upper panel) and the mass accretion rate (the lower 

panel). Due to this strong dependence on M, the energy deposition by the neutrino 

pair-annihilation would be important only for a phase in which L^tot ^ 10 51 ergs/s 

(see Figs. [TlT b) and fT2T b)). Figure fT2Tb) shows that for UNlOO-rigid, the duration 

of the neutrino emission (in the BH phase) with L^tot ^ 10 51 ergs/s is longer than 

1 s, and thus, a long-term energy deposition for a LGRB may be explained. Taking 

into account the dependence of the neutrino pair annihilation rate on the geometry 
of the torus,E3'E3KE3'E3 ^_ would be given hy m 

• 48 f lOOkm ^ f0.l\ 2 ( E v + E 9 \ 
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Fig. 14. Time evolution of the irreducible mass of the BH (the upper panel) and the mass accretion 
rate (the lower panel) for UNlOO-rigid. 
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where -Rf un and 0f un are the characteristic radius and the opening angle of the funnel 
region. denotes the collision angle of the neutrino pair. Thus a low-luminosity 
LGRB could be explained. 

In the HMNS phase, by contrast, the neutrino luminosity is huge as L u > 10 53 
ergs/s (see Figs. fTTT a) and fT2T a)). and hence, the deposition rate would be very large 
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If the outflows launched due to the mass accumulation mechanism can penetrate the 
stellar envelope, a system composed of a long-lived HMNS and a geometrically thick 
torus may be a promising candidate of the central engine of LGRBs a relatively short 
duration. 
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§4. Binary neutron star merger 

Coalescence of binary neutron stars (BNS) is one of the most promising sources 
for next-generation kilo-meter-size gravitational-wave detectors J3S>ED and also a 
possible candidate for the progenitor of SGRBsP^'M 1 Motivated by these facts, 
numerical simulations have been extensively performed for the merger of BNS in the 
framework of full general relativity in the past decade since the first success 8 -^* in 
2000 (see also, e.g., Refs. 33) for reviews). 

BNS evolve due to the gravitational radiation reaction and eventually merge. 
Before the merger sets in, each neutron star is cold (i.e., thermal energy of constituent 
nucleons is much smaller than their Fermi energy) , because the thermal energy inside 
the neutron stars is significantly reduced by neutrino and photon emissions^ in the 
long-term inspiral phase (typically > 10 8 yeard®) until the merger. By contrast, 
after the merger sets in, shocks are generated by hydrodynamic interactions. In 
particular, when a HMNS is formed in the merger, spiral arms are developed in its 
envelope and continuous heating occurs due to the collision between the HMNS and 
spiral arms (e.g., Refs. 184")) .185 ]) .186 ]) ). Newtonian simulations indeed suggest that 
by this process the maximum temperature increases to ~ 30-50 MeV, and hence, 
copious neutrinos are emittedPD'^'^''® Thus, to accurately explore the merger 
process, the evolution of the hot HMNS, and possible subsequent formation of a 
BH with a physical modeling, numerical-relativity simulations have to be performed 
incorporating microphysical processes such as neutrino emission and equation of 
state (EOS) based on a theory for the high-density and high-temperature nuclear 
matter. However, such simulations have not been done in full general relativity until 
quite recently (but see Ref. [19|) for a work in an approximate general relativistic 
gravity with finite-temperature EOS). Incorporation of microphysical processes is in 
particular important for exploring the merger hypothesis of SGRB because it may 
be driven through pair annihilation of neutrino-antineutrinos pairsP^'' 8 -^' 

In this section, we review our first results of numerical-relativity simulations for 
the BNS merger presented in Refs. l27|)J28p . which are performed incorporating both a 
finite-temperature EOSP'^ and neutrino cooling.^ In the following, we summarize 
the possible outcome formed after the merger, criteria for the formation of HMNS 
and BH, thermal properties of the HMNS and torus surrounding the formed BH, 
and neutrino luminosity and gravitational waveforms from the HMNS and in the 
BH formation. 

4.1. Initial condition and grid setting 

In Refs. I2"7 |) . l28|) . we focused only on the merger of equal- mass BNS, because 
the mass difference for the observed BNS is not very 

largeP^SJ To dat6) 

we have 

performed simulations for 5 models: For Shen-EOS, we employed three masses for 
each neutron star: Mns = 1-35, 1.5, and 1.6Mq (Mns is the gravitational mass of 
a neutron star in isolation). We refer to each model as models S135, S15, and S16, 
respectively. For Hyp-EOS, we employed Mns = 1-35 and 1.5M Q , and refer to two 
models as H135 and H15, respectively. The simulations were performed with the 
initial condition of about 3-4 orbits before the onset of the merger, until the system 
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relaxes to a quasi-stationary state. Quasi-equilibrium states of BNS were prepared 
as the initial conditions, as in Refs. [84"]).l85p. using the LORENE libraryPD 

There are two possible fates® of BNS: If its total mass M is larger than a 
critical mass M c , a BH will be formed soon after the onset of the merger, while a 
differentially rotating HMNS will be formed for M < M c . The value of M c depends 
strongly on the EOS. Because Shen-EOS is quite stiff, M c is much larger than the 
typical total mass of BNS, ~ 2.7 M Q ^®® as shown in Ref. EZD and below. Thus, 
with this EOS, the HMNS will be the frequent outcomes, as in the cases of stiff EOS 
with which M max > 2MqW^ By contrast, Hyp-EOS is not stiff in particular for a 
high-density range. Thus, a BH is often formed with this EOS, although a HMNS 
could be a transient outcome soon after the onset of the mergerPHJ 

Numerical simulations were performed preparing a non-uniform grid as in Ref. [85]). 
The inner domain was composed of a finer uniform grid and the outer domain of a 
coarser nonuniform grid. The grid resolution in the inner zone is chosen so that the 
major diameter of each neutron star in the inspiral orbit was covered by 60 and 80 
grid points for low- and high-resolution runs, respectively: We always performed sim- 
ulations for both grid resolutions to confirm that the convergence, sufficient to draw 
a scientific conclusion on the final outcome, gravitational waveforms, and neutrino 
luminosities, is approximately achieved. Outer boundaries are located in a local wave 
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Fig. 16. The rest mass of a torus surrounding the formed BH as a function of time for models H135, 
H15, and S16. £bh denotes the time at the formation of the BH. 

zone (at ~ 560-600 km along each coordinate axis which is longer than gravitational 
wavelength in the inspiral phase). During the simulations, we checked the conserva- 
tion of the baryon rest-mass, total gravitational mass (Arnowitt-Deser-Misner mass 
plus radiated energy of gravitational waves), and total angular momentum (includ- 
ing that radiated by gravitational waves), and found that the errors are within 0.5%, 
1%, and 3%, respectively, for the high-resolution runs within the physical duration 
~ 30 ms. 

4.2. Merger and subsequent evolution 

Figure [TBI plots the maximum rest-mass density, p ms , x , maximum matter temper- 
ature, T max , and maximum hyperon fraction in mass Xj\ <max as functions of t — t meTge 
where t m erge is the onset time of the merger. For t < t merge , /3 max is approximately 
constant besides a small decline due to tidal elongation, while for t > i merge , it 
gradually increases because a HMNS is formed at least temporarily irrespective of 
models, and subsequently contracts due to the dissipation of the angular momentum 
by the gravitational-wave emission. Thus, /5 max increases in the gravitational radi- 
ation time scale. The subsequent evolution process depends on the mass and EOS. 
For models S135 and S15, the degree of non-axial symmetry of the HMNS becomes 
low enough at t — t mcrgc ~ 20 ms that the emissivity of gravitational waves is sig- 
nificantly reduced. Because no dissipation process except for the neutrino cooling is 
present, the HMNS will be alive at least for the cooling time scale before collapsing 
to a BH (see below). For models S16, H135, and H15, the HMNS collapse to a BH 
at t — tmerge ^ 10 ms after the gradual contraction due to the gravitational-wave 
emission and a massive disk of ~ 0.03-0.1-Mq is formed around the BH. It should be 
noted that for Hyp-EOS, a BH is formed at t — t meTge ~ 10 ms even with the total 
mass 2.7 'M@. This is due to the softening effect by the appearance of A hyperons: 




36 Y. Sekiguchi, K. Kiuchi, K. Kyutoku and M. Shibata 




-3£3 -2E -113 □ ID 33 3D -3G -213 -113 □ 1E3 33 3D 



X[km] 

Fig. 17. Contour maps in the x-z plane of the rest-mass density (top left), the electron fraction 
(top right), the temperature (bottom left), and the total neutrino emissivity (bottom right) at 
t 16.7 ms after the onset of the merger for model S135. 

See the bottom panel of Fig. [15j which shows that XA,max. increases steeply just 
before the HMNS collapses to a BH. 

The evolution of T max plotted in Fig. 1151 shows that the HMNS formed just after 
the merger are hot with T max ~ 50-70 MeV (much higher than that in the ordinary 
supernova and as high as that in the HMNS formed after the collapse of the massive 
stellar core; cf. §[3|). Such a high temperature is achieved due to the liberation of 
the kinetic energy of the orbital motion at the collision of two neutron stars. For 
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Fig. 18. The same figure as Fig. [17] but in the x-y plane. 



the case that a long-lived HMNS is formed, subsequently, T max decreases due to the 
neutrino cooling, with the maximum luminosity 3-10 x 10 53 ergs/s (see Fig. [23]), but 
relaxes to a high value with 25-50 MeV when the HMNS relaxes to a quasi-steady 
state. Around the HMNS, spiral arms are formed and shock heating continuously 
occurs when the spiral arms hit the HMNS (see Figs. [T7H22I for snapshots). Due to 
this process and because of the long neutrino-cooling time scale, the temperature 
(and thermal energy) does not significantly decrease in ~ 100 ms: We estimated the 
cooling time scale as E t h/L u ~ 2-3 s where E^ is the total thermal energy of the 
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Fig. 19. Contour maps of the rest-mass density for a HMNS phase at t ~ 17.5 ms after the 
merger (left panels) and that for a BH phase at t « 26.8 ms after the onset of the merger for 
model H135. The upper and lower panels show the configuration in the x-y and x-z planes, 
respectively. The blue circle of the right panels shows the location of the apparent horizon. 



HMNS. 

For the case that a BH is eventually formed, the maximum temperature raises 
significantly to > 100 MeV just before the BH formation. This is simply due to the 
adiabatic compression effect. For models S16, H135, and H15, a torus surrounding 
the BH is subsequently formed. The typical maximum density and temperature of 
the torus are ~ 10 13 g/cm 3 and 20 MeV, respectively, with the mass ~ 0.03-0. 1M@ 
(see Fig. [To]) . This mass has a correlation with the lifetime of the HMNS; for the 
longer lifetime (e.g., for model S16), the torus mass is larger (compare the mass of 
the torus in Fig. [To]) . The reason is that during the evolution of the HMNS which 
is deformed in a non-axisymmetric manner, the angular momentum is transported 
from the inner to the outer region via the hydrodynamic torque associated with 
its non-axisymmetric structure. Thus, the longer lifetime helps increasing the mass 
element of a sufficiently large specific angular momentum which can escape falling 
into the formed BH. This fact implies that stiff EOS are favored for the formation 



Current Status of NR simulation 



39 




X [km] 

Fig. 20. The same figure as Fig. [19] but for the electron fraction, 
of a massive torus. 

Figures [T7] and [18] plot the contour maps of the rest-mass density, electron frac- 
tion, matter temperature, and total neutrino luminosity of a HMNS for model S135 
ms in the x-z and x-y planes, respectively, at which it already 
relaxed to a semi-final quasi-steady state. This shows that the HMNS is weakly 
spheroidal and the temperature is high (T ~ 30 MeV) in its outer region. The 
neutrino luminosity is also high in its outer region, in particular, near the polar sur- 
face. With the fact that the rest-mass density is relatively small near the rotation 
axis above the polar surface, this is a favorable feature for the merger hypothesis of 
SGRB; pair annihilation of neutrinos and anti-neutrinos could supply a large amount 
of thermal energy which may drive a fire ball along the rotation axis. As in § 13.4] 
the neutrino pair annihilation rate is estimated as 

■ ln51 . / , !01mW0.3\ 2 /E„ + £A 
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Fig. 21. The same figure as Fig. [19] but for the temperature. 

which would be sufficient for driving SGRBs. The pair annihilation efficiency has 
been approximately estimated in the previous works^'^'^'^ and our result is 
consistent with these works. 

Possible reasons that HMNS are formed are; (i) it is rapidly rotating with the 
period ~ 1 ms, and hence, the centrifugal force increases the possible mass that can 
be sustained; (ii) because it is hot, the thermal energy enhances the pressure. We 
find that the rotational velocity with the period ~ 1 ms does not play a substantial 
role. Exploring in detail Shen-EOS for the high density tells us that the effect of 
the thermal energy is significant and can increase M max by ~ 20-30% for a high- 
temperature state with T > 20 MeV. This indicates that the HMNS will alive before 
collapsing to a BH for a long cooling time > 1 s. At the time when the HMNS collapse 
to a BH, it will be close to a spherical configuration with low temperature due to 
long-term gravitational-wave and neutrino emissions. Thus, observable signals from 
the late-time collapse will not be remarkable. 

Figure [23] plots neutrino luminosities as functions of time for three flavors (u e , P e , 
and sum of v x ). It is found that electron anti-neutrinos are dominantly emitted for 
any model. The reason for this is as follows: The HMNS has a high temperature, and 
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Fig. 22. The same figure as Fig. [19] but for the total neutrino emissivity. 



hence, electron-positron pairs are efficiently produced from thermal photons, in par- 
ticular in its envelope. Neutrons efficiently capture positrons to emit anti-neutrinos 
whereas electrons are not captured by protons as frequently as positrons because 
the proton fraction is much smaller. This hierarchy in the neutrino luminosities was 
reported also in Refs. l2T ]hl22"l) . 

Soon after the BH formation for models S16, H135, and H15, /jl/t neutrino 
luminosity steeply decreases because high-temperature regions are swallowed into 
the BH, while luminosities of electron neutrinos and anti-neutrinos decrease only 
gradually because these neutrinos are emitted via charged-current processes from the 
massive accretion disk. We here note that magnetic fields, which are not taken into 
account in the present simulations, could be amplified significantly in the accretion 
disk 8 ^ 1 and may play a role in the late evolution of the BH-disk system. 

The anti-neutrino luminosity for the long-lived HMNS is Lp ~ 1.5-3 x 10 53 ergs/s 
with a small time variability. It is by a factor of ~ 1-5 larger than that from proto- 
neutron stars formed after supernovaeP while it is comparable to the luminosities 
for UNlOO-rigid model. The averaged neutrino energy is ep ~ 20-30 MeV. The 
sensitivity of water-Cherenkov neutrino detectors such as Super-Kamiokande and 
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future Hyper-Kamiokande have a good sensitivity for such high-energy neutrinos 
in particular for electron anti-neutrinos.^ The detection number for electron anti- 
neutrinos is approximately estimated by a AT Ly / {AttD 2 e^) where a is the total cross 
section of the detector against target neutrinos, AT is the lifetime of the HMNS, and 
D is the distance to the HMNS. For a one-Mton detector such as Hyper-Kamiokande, 
the expected detection number is > 10 for D < 5 Mpc with AT ~ 2-3 s, based on an 
analysis of Ref . [93]) . Thus, if the BNS merger fortunately happens within D ~ 5 Mpc, 
neutrinos from the HMNS may be detected and its formation may be confirmed. Note 
that gravitational waves from the HMNS will be simultaneously detected for such a 
close event (see below), reinforcing the confirmation of the HMNS formation. 

4.3. Gravitational waves 

Figure 1247 a) plots the plus mode of gravitational waves as a function of 
*ret — Emerge where t Tet is the retarded time, t ret = t — D — 2Mlog(D/M) (M = 
2Mns)- Gravitational waves are extracted from the metric through the outgoing 
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Fig. 24. (a) Gravitational waves observed along the axis perpendicular to the orbital plane for the 
hypothetical distance to the source D = 100 Mpc for all the models, (b) The effective amplitude 
of gravitational waves as a function of frequency for D = 100 Mpc. The noise amplitudes of a 
broadband configuration of Advanced LIGO (bro. LIGO) and KAGRA are shown together. 
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Fig. 25. /gw(*) in the HMNS evolution phase, smoothed by a weighted spline, for models H135, 
S135, and S16. 



component of the complex Weyl scalar, ^4, in the local wave zone. The waveforms 
are composed of the so-called chirp waveform, which is emitted when the BNS is in 
an inspiral motion (for t ret < Emerge ), and the merger waveform (for t ret > i mer ge)- 
Gravitational waves from the inspiral phase (for t ret < Emerge) agree well with each 
other for the models with Hyp-EOS and Shen-EOS for the same mass. On the 
other hand, quasi-periodic gravitational waves from the HMNS (for t ie t > Emerge) 
show several differences. First, the amplitude of quasi-periodic gravitational waves 
damps steeply at the BH formation for H135 and H15. This is because the HMNS 
collapse to a BH before relaxing to a stationary spheroid. Second, the characteristic 
gravitational- wave frequency, /gWj increases with time for Hyp-EOS models, while 
it is approximately constant for Shen-EOS models with / pea k ~ 2.0-2.5 kHz which 
depends weakly on M. These facts are clearly observed in the effective amplitude 
(see Fig. 121Tb)) defined by h e s{f) = 0Af\h(f)\ where h(f) is the Fourier transform 
of h + — ih x with h x being the cross mode and the factor 0.4 comes from taking the 
average in terms of the random direction to the source and rotational axis of the 
HMNS. Reflecting a shorter lifetime of the HMNS in Hyp-EOS models, the peak 
amplitude of /i e fr(/) is smaller, in particular for H15 where the HMNS survives only 
for a short period ~ 3 ms. Reflecting the shift of the characteristic frequency, the 
prominent peak in /i e g for Hyp-EOS models (H135 and H15) is broadened. The 
reason for this is described as follows in more detail. 

In the case that hyperons are absent, the HMNS slightly contract during their 
evolution simply due to the angular momentum loss (weakening centrifugal force). 
By contrast, in the case that hyperons are present, X\ increases with the contraction 
of the HMNS, resulting in the relative reduction of the pressure. As a result, the 
HMNS contracts by a larger fraction. Recent studies showed that /gw is associated 
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with the frequency of an /-mode which is approximately proportional to yMn/i?g 

where Mh and Rn are the mass and radius of the HMNSPl This indicates that 
/gw should increase with time. To see that this is indeed the case, we show /gw(= 
d(p NP /dt) calculated from W A = |Psi 4 |e i<?w in the HMNS phase for H135, S135, and 
S16 in Fig. 1251 It is clearly seen that the mean value of /gw is approximately constant 
for Shen-EOS models; /gw ~ 2.1 and 2.5 kHz for S135 and S16, respectively. By 
contrast, /gw for H135 increases with time (from /gw ~ 2.0 kHz at t Tet — i mer ge = 2 
ms to ~ 2.5 kHz at t re t — tmerge = 10 ms) as the HMNS becomes compact. 

Figure 124( b) shows that if a HMNS with the lifetime > 10 ms is formed, the 
effective amplitude h e g ~ 4-6 x 10~ 22 at D = 100 Mpc for Shen-EOS models and 
2 x 10~ 22 at D = 100 Mpc for H135 with / peak « 2.0-2.5 kHz which depends weakly 
on M and EOS. For the detection, Hyp-EOS is obviously unfavored. By contrast, 
for Shen-EOS, the maximum amplitude for a hypothetical distance of 100 Mpc is as 
high as the sensitivity curve of a specially-designed version of advanced gravitational- 
wave detectors such as broadband LIGO,'^' which has a good sensitivity for a high- 
frequency band. This suggests that gravitational waves from the HMNS oscillations 
may be detected with S/N = 5 if D < 20 Mpc. If the source is located in an 
optimistic direction, the detection with S/N = 5 may be possible for D = 50 Mpc. 

§5. Summary 

We have described our latest results of numerical-relativity simulations of ro- 
tating stellar core collapses to a BH and BNS mergers, performed incorporating a 
finite-temperature EOS and neutrino cooling effects. The following is the summary 
of our latest findings and prospects for the near future. 

5.1. Stellar core collapse 

We presented our latest results of axisymmetric simulations of very massive stel- 
lar core collapsing to a system composed of a rotating BH and surrounding disk/torus 
in full general relativity. The simulation were performed taking into account of the 
microphysical processes and the neutrino cooling. Because progenitor models of 
LGRBs suggested in the literatures®-' propose a possibility that they may have an 
entropy higher than that of ordinary supernova cores, we employed a model of a pre- 
supernova core with a high entropy of s/fcg =4 calculated by Umeda and Nomoto®' 1 
together with hypothetical two rotational profiles (UNlOO-rigid and UNlOO-diff). 

As in the collapse of ordinary supernova cores, the gravitational collapse sets 
in due to the photo-dissociation of heavy nuclei and the electron capture. The 
collapsing core eventually experiences a core bounce, forming a shock wave. Then 
a HMNS, which is supported by the centrifugal force and the thermal pressure is, 
formed. The neutrino luminosity in this HMNS phase is larger than than that in 
the ordinary supernova as > 5 x 10 53 ergs/s. The HMNS eventually collapses to 
a BH irrespective of the initial rotational profiles. However, the dynamics and the 
geometry of the final outcome depend strongly on the degree of the initial rotation. 

For the model UNlOO-rigid, the shock wave formed at the core bounce is de- 
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formed to be a torus-like shape. Then the infalling materials are accumulated in the 
central region after they pass through the oblique shock formed at the tours-shaped 
outer region of the HMNS. As a result, the thermal energy is efficiently stored at the 
surface of the HMNS due to the dissipation of the kinetic energy of the accumulated 
materials, driving outflows. After the collapse of the HMNS, a torus is formed around 
the BH. We found that the torus shows a time variability. The total neutrino lumi- 
nosity emitted from the torus around the BH amounts to L v t Q t ~ 

10 5i_ 1Q 52 ergs / S) 

which lasts for a long duration of > 1 s. Associated with the time variability of 
the BH-torus system, the neutrino luminosities also show a violent time variability. 
Such a long-term high luminosity with the time variation may be related to the time 
variability that LGRBs show. For the model UNlOO-diff, by contrast, a geometri- 
cally thin disk is formed around the BH and the BH-disk system shows essentially no 
time variability. Remarkably, the above differences in the dynamics and the outcome 
stem from a small difference in the initial rotational profile. 

We also calculated the characteristic gravitational-wave strain /i c har for UN100- 
rigid. The effective amplitude is as large as ~ 10 -20 at / ~ 1 kHz for a hypothetical 
event occurred at a distance of 10 kpc. This shows a possibility that we may observe 
multi-messenger information, namely, gravitational waves, neutrinos, and electro- 
magnetic radiation from such a nearby event and may obtain a clue to understand 
the stellar core collapse and the central engine of LGRBs. 

5.2. Binary neutron star merger 

We showed that for a stiff, purely nucleonic EOS, a HMNS is the canonical 
outcome and a BH is not promptly formed after the onset of the merger as long as 
the total mass of the system is smaller than 3.2M . The primary reason is that 
the thermal pressure plays an important role for sustaining the HMNS. We further 
showed that the lifetime of the formed HMNS with mass < 3M© would be much 
longer than its dynamical time scale, i.e., » 10 ms, and will be determined by the 
time scale of the subsequent neutrino cooling. The neutrino luminosity in the early 
evolution phase of the HMNS was shown to be high as ~ 3-10 x 10 53 ergs/s. The 
effective amplitude of gravitational waves averaged over the random source direction 
and orbital plane inclination is h e g = 4-6 x 10 -22 at / pca k = 2.1-2.5 kHz for a 
hypothetical source distance of D = 100 Mpc. If the BNS merger happens at a 
relatively short source distance ~ 20 Mpc or is located in an optimistic direction with 
D ~ 50 Mpc, such gravitational waves may be detected by advanced gravitational- 
wave detectors with S/N = 5, and the HMNS formation will be confirmed. 

For an EOS in which effects of hyperons are taken into account, the EOS becomes 
softer than the purely nucleonic EOS. With this EOS, a BH is often formed in a short 
time scale after the onset of the merger, although a HMNS could be a transient 
outcome with a short lifetime < 10 ms. Because the EOS becomes soft during the 
evolution of the HMNS, the compactness significantly changes in a short time scale 
in this EOS. This is well reflected in gravitational waveforms and their spectra. 
Specifically, the characteristic frequency changes with time. This effect reduces the 
amplitude at a peak frequency of gravitational waves in the Fourier space, and make 
a feature unfavorable for the detection of gravitational waves. Roughly speaking, 
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the allowed distance for the detection of gravitational waves from the HMNS is by a 
factor of 2 smaller than that in the nucleonic EOS for the same mass of BNS. 

5.3. Future prospects 

5.3.1. Massive stellar core collapse 

As mentioned in § [3l we did not take into account effects of the neutrino heating 
in the numerical simulations to date. Recently, we have developed-® a formulation 
for numerical simulations of general relativistic radiation transfer based on Thome's 
moment formalism!^* Based on this formalism, we have already performed general 
relativistic radiation magnetohydrodynamics (GRRMHD) simulations^ f Q r the evo- 
lution of a system composed of a BH and a surrounding torus with a simplified treat- 
ment of microphysics, as a step toward a more physical modelling. Furthermore, we 
have succeeded in implementing a code which can solve the neutrino transfer with 
a detailed microphysics (in preparation). Using this code, we plan to perform simu- 
lations of the stellar core collapse to explore a supernova explosion mechanism and 
the formation of a BH in full general relativity. 

Throughout the study of the massive stellar collapse in this article, we assume 
the axial symmetry in the simulations. However, this assumption would be invalid 
if non-axisymmetric instabilities set in. For example, in Ref. [96]), it has been shown 
that the BH-torus systems could exhibit the so-called Papaloizou-Pringle instabil- 
ity.^ Also, for a rapidly rotating HMNS, non-axisymmetric instabilities such as 
a bar-mode instability could set i n .USD .lifloj ,[ioiji ,ETJ) Q nce these instabilities turn 
on, the torus/HMNS may deform to be a highly non-axisymmetric structure. This 
will enhance the angular momentum transport in the torus and HMNS, and the 
evolution processes of these systems may be modified. We plan to perform a three 
dimensional simulation to explore if the non-axisymmetric instabilities set in and 
play an important role for the collapsar models. 

5.3.2. Binary neutron star merger 

To date, we performed the BNS merger simulations for the case of the equal-mass 
binaries. As a straightforward extension of the previous studies, we plan to perform 
simulations for unequal-mass binaries, for which the merger dynamics, gravitational 
waveforms, and the mass of the disk will be modified. We also plan to perform merger 
simulations of BHNS binaries. It has been reported (see Ref. HU2p and references 
therein) that a massive disk of Afdisk O.IMq can be formed for stiff EOSs even 
when the BH rotates moderately (obh ^ 0.5). Such a system is a promising candidate 
of the central engine of SGRBs. Only simulations with detailed microphysics will 
enable a quantitative study for this merger hypothesis of SGRB. 

In the study of the BNS merger simulation, we totally ignore the effect of mag- 
netic fields. Recent studies on magnetized BNS merger showed that the magnetic 
fields in the merger and post-merger phases have an impact on the dynamics of the 
torus formed around the BHP^ This is because the angular velocity inside the torus 
has a steep gradient. Thus, the magnetic field is subject to the amplification via 
magnetic winding and/or magneto-rotational instability This condition holds in 
the HMNS because it has strong and rapid differential rotation as discussed in § 4.2. 
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We plan to incorporate magneto-hydrodynamics in our code and explore its impact 
on the evolution of HMNSs. 

Acknowledgments 

We thank H. Umeda for providing us the presupernova model (UN100) adopted 
in this work. Numerical simulations were performed on SR16000 at YITP of Kyoto 
University, on SX9 and XT4 at CfCA of NAOJ, and on the NEC SX-8 at RCNP in 
Osaka University. This work was supported by Grant-in-Aid for Scientific Research 
(21018008, 21105511, 21340051, 21684014, 22740178, 23740160), Grant-in-Aid on 
Innovative Area (20105004), and HPCI Strategic Program of Japanese MEXT. 

Appendix A 

Electron and positron captures 

In this section, we briefly summarize our methods of handling electron and 
positron captures based on Ref . 15?]) . and give the explicit forms of 7^, 7p e °, Q^ e , and 
Ql c e in Eqs. d755|> . (12^MD . (I2~^dD . and (l^oTD . for completeness 

A.l. The electron and positron capture rates 7^° and 7p e c 

The 'net' electron fraction is written as Y e = Y_ — Y + where Y— (Y + ) denotes 
the number of electrons (positrons) per baryon including pair electrons. Then the 
electron-neutrino number emission rate by the electron capture and the electron- 
anti-neutrino number emission rate by the positron capture are given by 

7 iocai = _y_ = _ { y[ + yh), (A-l) 

7 iocai = _ y+ = _ ( y/ + y^ (A . 2) 

where the electron and positron capture rates are decomposed into two parts, capture 
on by free nucleons (with the superscript /) and on heavy nuclei (with the superscript 
h). In the following, we will present the explicit forms of Y_, Y^, Yh, and Y+. 

A. 2. Capture on free nucleons Y* 

The electron capture rate (including the contribution of the inverse reaction of 
the neutrino capture) on free nucleons (y£) is given by 

YL = X n \ UeC 'f - X p X ec ' f , (A-3) 

where A ec '^ is the specific electron capture rate on free protons, \ UeC >f is the specific 
electron-neutrino capture rate on free neutrons, and X p and X n are the mass fraction 
of free protons and neutrons, respectively. Based on a balance argumentp' 1 one can 
show that \ UeC >f is related to A ec '^ by 

^ eCj = exp U e ~Ve- ^) X ccJ , (A-4) 

where r\ Ve and r/ e are the chemical potentials of electron neutrinos and electrons in 
units of JzbT and 5m = (m n — m p )c 2 . Furthermore, we use the following Saha's 
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relation for non-degenerate free nucleons, 

(Sm \ 
T]n ~ rip+ k^f) ' ( ~ A ' 5 ' ) 

where r/ n and ij p are the chemical potentials of free neutrons and protons in units of 
kpT. Then we obtain 

Yl = [exp (rj Ve -n. + rjn- r, p ) - 1} X p X cc ^ '. (A-6) 

The positron capture rate (including the contribution of the inverse reaction) on 
free nucleons is similarly given by 

Y( = X p X p ^f - X n \^f = [exp ( VPe + Ve + ri p - %) - 1] X n \^ , (A-7) 

where r]p e is the chemical potential of electron-anti- neutrinos in units of kpT, A pc is 
the specific positron capture rate on free neutrons, and \ VeC ^ is the specific electron- 
anti- neutrino capture rate on free protons. 

A. 3. Capture on heavy nuclei Y h 

The electron capture rate (including the contribution of the inverse reaction of 
the neutrino capture) on a heavy nucleus of mass number A (Y^) is given by^ 

yh = ^D_ x u e c,h _ ^£\ ec ' h , (A-8) 

where \ ec ' h is the specific electron capture rate on the parent nucleus (mass fraction 
Xp), \ UeC > h is the specific electron- neutrino capture rate on the daughter nucleus 
(mass fraction Xd), and A is the atomic mass of the parent and daughter nuclei. 
In the present simulations, we set Xp> = Xp = Xa- Then, under the assumption 
of a nuclear statistical equilibrium, one may approximate the capture rate on heavy 
nuclei asJ23 

Y h « [exp (r, Ue - Ve + r, n - Vp ) - 1] ^A\^\ (A-9) 

Similarly, the positron capture rate (including the contribution of the inverse 
reaction) on heavy nuclei (Yi) is given by 

yh = XDyv^h _ Xp_ xV c,h „ [exp ^ +r]e+7]p _ %) _ !] ^APC/*. (A . 10) 

A. 4. The specific capture rate A 

The specific electron and positron capture rates on free nucleons and on heavy 
nuclei are written in the same form asr^ 

A ec,/ = J^_jocJ A pcJ = J^_jpcJ (A . n) 

X ec,h = ln2 jec,h X pc,h = ln2 pc,h (A . 12) 
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where I ec, f and J pc '^ are the phase space factors for the electron and positron cap- 
tures on free electrons, and I ec,h and iv c > h are those on heavy nuclei. (/i) e ff's are 
the effective ft- values introduced by Fuller et al.jS3 which is essentially the same as 
the square of the nuclear transition matrix. 
The phase space factors are given by 

irfi ■ C vj 'r- • — i 
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dr], (A- 13) 

dr), (A-14) 

drj, (A-15) 

dr], (A-16) 



where C ec '^! C pc '^> C ec,h i an d ( pc,h are the nuclear mass-energy differences for the 
electron and positron captures in units of k B T. The superscripts 'f and 'h' again 
denote free nucleons and heavy nuclei. The nuclear mass-energy differences for the 
capture on free nuclei are given by 

C' f = -CF' f *V P -Vn. (A-17) 

We follow Fuller et alP3 for the nuclear mass-energy differences in the capture on 
heavy nuclei: In the case of N < 40 or Z > 20 (referred to as 'unblocked' case), we 
set 

C e C ,h = _ ec ,h^ r]p _ rjn (A . 18) 

In the case of N > 40 or Z < 20 (referred to as 'blocked' case), on the other hand, 
we set 

5MeV 

Vp ~ Vn 



c 



pc,/i 



-Vp + Vn + 



k B T ' 
5MeV 
k B T ' 



(A49) 
(A-20) 



Then, the threshold value of the electron and positron captures is given by r/o = 
m e c 2 /(k,BT) for £ > —m e c 2 /(kBT) and tjq = |£| for ( < —m e c 2 /(kBT) where we 
have dropped the superscripts 'ec', 'pc', '/', and '/i' in £ for simplicity. 

The effective ft- value of the electron or positron capture on free nuclei is given 
by (e.g. Ref. EZD 

logical = logio(/i) e P ff' / ~ 3.035. (A-21) 

We follow Fuller et alW^ for the effective /t-value of the capture on heavy nuclei, 
who proposed to use 

unblocked r] e < \( ec > h \ 
unblocked r, e > \( ec ' h \ , (A- 22) 
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log l0 {ft)^ h 
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3.2 unblocked rj e < |C pc,/l | 
2.6 unblocked rj e > |C pc ' h | , 
2.6 + ffp blocked 



(A-23) 



where T 9 = T/^IO 9 ^). In this expression, the thermal unblocking effect?® is readily 
taken into account. In the thermal unblocking, it costs ~ 5.13 MeV to remove a 
neutron from a filled orbital I/5/2 and place it in the gd-shellW^ 

A. 5. Energy emission rates and Qp^ 

The neutrino energy emission rates associated with the electron and positron 
captures in units of m e c 2 s _1 are given by^ 
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dr/, (A-25) 
dr/. (A-26) 



In Eqs. (|A-24p - (|A-26p . we have dropped the superscripts '/' an d 'h' 1 in 7r ec , 7r pc , J ec , 
d pc , (ft)th (/*>eff> C ec , and C pc for simplicity. 

The average energy of the electron neutrinos produced by electron and positron 
captures is defined, in units of m e c 2 , as 
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Then, the local neutrino energy emission rates by the electron and positron captures 
per unit volume is given by 
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Appendix B 

- Neutrino pair processes 



(A-28) 
(A-29) 



In this section, we briefly summarize our method of handling the pair processes 
of the neutrino emission and give the explicit forms of 7 P ^, l^ e , 7^, r p e mS ! 7 pa p > 
ll% 7?X S , Ql% Ql% Qu T JT> Ql%i QtL ^d for completeness. " 

B.l. Electron-positron pair annihilation 

We follow Cooperstein et alPO 1 for computing the rate of neutrino emission by 
the electron-positron pair annihilation. The number emission rate of v e or v e by the 
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electron-positron pair annihilation can be written as 

iriu C^p] a c (k B T) 



Pair 



F 3 (r7 e )F 3 (-r ?e )(block)PJ e , 



(B-l) 



p 36vr 4 m e ¥ (He) 6 

where a ~ 1.705 x lO^cm" 2 and C* % = (Cy - C A f + (Cy + C A ) 2 with Cy = 
g + 2sin 2 #iy and Ca = \- The Weinberg angle is given by sin 2 d\y « 0.23. Using 
the average energy of neutrinos produced by the pair annihilation, 
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the blocking factor (block) {^f is evaluated as 
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The associated neutrino energy emission rate by the pair annihilation is given by 

P 
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(B-4) 



Similarly, the number emission rate of v x or v x by the electron-positron pair 
annihilation and the associated energy emission rate are given by 
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where C l/x p x = (Cy — Ca) 2 + (Cy + d — 2) 2 . The average neutrino energy and the 
blocking factor are given by 
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where r\ Vx = r) Vx because they are produced only by the pair processes. 

B.2. Plasmon decay 

We follow Ruffert et alPD for computing the pair creation rate of neutrinos by 
the decay of transversal plasmons. The number emission rate of v e or v e can be 



written as 
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where a^ nc 1/137 is the fine-structure constant and j p « 2y/ (afi ne /97r)(7r 2 + 3rj e ). 
The blocking factor is approximately given by 
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is the average energy of neutrinos produced by the plasmon decay. The associated 
neutrino energy emission rate is given by 



^plas _ P plas / _ v plas 
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Similarly, the number emission rate of v x or v x by the plasmon decay and the 
associated energy emission rate are given by 

■US. = -fc^S^T^d + 7 r )(bloc k ,S„ (B-13) 
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where the average neutrino energy is (e^p^) 131 ^ = {(-v e v e ) pXas and the blocking factor 
is given by 
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B.3. Nucleon-nucleon bremsstrahlung 

We follow Burrows et alP^ for computing the pair creation rate of neutrinos 
by the nucleon-nucleon bremsstrahlung radiation. They derived the neutrino energy 
emission rate associated with the pair creation of v x or v x by the nucleon-nucleon 
bremsstrahlung radiation without the blocking factor as 
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where £ Brems ~ 0.5 is a correction factor and the average energy is 
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To obtain the 'blocked' neutrino energy emission rate we multiply the blocking factor, 



(block)?? 



Brcrns 

'x Vx 



1 + exp n, 



Brcms 



k B T 



-1 






1 + exp ^ 



( e VxVx) 



Brcms 



k B T 



(B48) 



54 



Y. Sekiguchi, K. Kiuchi, K. Kyutoku and M. Shibata 



to give 

Qu T JT = ^r U (block)^ 8 - 
The number emission rate of v x or f/^ is readily given by 

4.5 



7^ = 3.62 x 10 5 C Brems [Xl + Xl + fx n X p j m u p J (block)*— 

(B-20) 

Noting that the weak interaction coefficients of the bremsstrahhmg radiation 
are^3 (1 - Cy) 2 + (1 - Ca) 2 for the pair creation of v x v x and C v + C\ for the pair 
creation of v e v e , the number emission rate and the associated energy emission rate 
for v e or v e are written as 
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Appendix C 

Neutrino diffusion rates 

We follow Ref. |22|) for computing the diffusive neutrino-number emission rate 
7^? and the associated energy emission rate Q(y, in Eqs. (|2-54p and (|2-55p . An 
alternative definition of the diffusion rates are found in Ref. I2T]) . 

C.l. Neutrino diffusion rates 

To calculate the neutrino diffusion rates 7^? and Qjyj > we first define the neu- 
trino diffusion time. In this paper, we consider cross sections for scattering on nuclei 
(o~l c A ), and on free nucleons (a^ p and <r^), as well as that for absorption on free 
nucleons (cr* and cr^). 

Ignoring the higher-order correction terms in neutrino energy E v , these neutrino 
cross sections can be written in general as 

a(E u ) = Eld, (C-l) 

where a is a 'cross section' in which E 2 dependence is factored out. In practice, the 
cross sections contain the correction terms which cannot be expressed in the form 
of Eq. (jC-ip . We take account of these correction terms, approximating neutrino- 
energy dependence on temperature according to 

E v « k B T^\. (C-2) 

The opacity is written as 
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and the corresponding optical depth is calculated by 

T (E V ) = j k{E v )qIs = E l j R ds = E 2 v t. (C-4) 
Then, we define the neutrino diffusion time by 

T^(E U ) = a^^^r(E v ) = E^t = E^, (C-5) 

C CK 

where the distance parameter Ax{E v ) is given by 

Ax{E u ) = ^|4_ (C-6) 

Note that T dlff can be calculated only using matter quantities. Here, a diS is a 
parameter which controls how many neutrinos diffuse outward and we chose it to be 
3 following Ref. |2T|) . For a larger value of a dlff , the corresponding neutrino emission 
rate due to diffusion becomes smaller. 

Finally, we define the neutrino diffusion rates by 

diff_™« f MEu) , F ^ 1 rn u 4TTcg u k 

difr _ /• E v n v {E v ) 1 47rcg„ k 2 

Q ^ = J T**(E V ) dEv ~ ^lhcfT* T FlM - (C ' 8) 

C.2. Summary of cross sections 

In this subsection, we briefly summarize the cross sections adopted in the present 
neutrino leakage scheme. 

C.2.1. Neutrino nucleon scattering 

The total v-p scattering cross section o~ p for all neutrino species is given by 

<p = T (^) 2 ^ Cv ~ 1 ) 2 +^\{Ca - I) 2 ] W; c (E u ), (C-9) 

where gA is the axial- vector coupling constant gA ~ —1.26. W^ c is the correction 
for the proton recoil. We use the exact expression derived by Horowit^® for high 
neutrino energies E v /m p c 2 > 0.01. However, the exact expression has a behavior 
which is inconvenient to treat numerically (such as 0/0). Thus we adopt expanded 
forms of the exact expression in low neutrino energies K/m p c 2 < 0.01, which give 

W1 C (E U )^ 1 - 1.524-^ + 1.451 (-%) , (C-10) 

y m p c z \ m p c z J 

W; c {E D )^l - 6.874-^ + 29.54 , (Oil) 



for neutrinos and anti-neutrinos, respectively. Note that in the case of black hole 
formation, the neutrino energy becomes large and this correction becomes important. 
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On the other hand, the total v — n scattering cross section a n is 

r) 2 [l + 3 5 i]^ c . (C-12) 



sc 

~ 16 



/TP \ 2 



1 2 



We evaluate W^ c by the same method as for Wp°. The expanded forms in the low 
neutrino energies are 

W1 C (E V ) « 1 - 0.7659-^ - 1.3947 ( , (C-13) 

m n c 2 \vn n & ) 

W s n c {E v )^ 1 - 7.366-^ +33.25 ( , (C-14) 

for neutrinos and anti- neutrinos, respectively. 
C.2.2. Coherent scattering of neutrinos on nuclei 

The differential cross section for the v-A neutral current scattering is written 



as 59 ' 



da B J (To / E u \ 2 2 o 



^ 2 [WC FF + C LOS ] (5io„)(l + cos#), (C-15) 



<ii7 647T \m e c 2 
where is the azimuthal angle of the scattering and 

W = l-— (l-2sin 2 ^)- (C-16) 

(<5ion)) Cl qs, and Cff are correction factors d ue t o the Coulomb interaction among 
the nucleij 1 ^* due to the electron polarization, 108 ' and due to the finite size of heavy 
nuclei Because it is known that the correction factor Clos is important only for 
low-energy neutrinos we consider only («Si on ) and Cff- 

The correction factor due to the Coulomb interaction among the nuclei is given 

by 

(S ion } = ^- J dcos9{l + cos9)(l-cos6)S ion . (C47) 
4 J-i 

Itoh et alP^S presented a detailed fitting formula for the correction factor. However, 
the fitting formula is so complicated that we use a simple approximation based on 
Ref-ECED, in which 

5ion w WT^hmf' (C ' 18) 

where q = (2E V / He) sin(0/2), ai = (47rn y i/3)~ 1 / 3 is the ion-sphere radius, ua is the 
number density of a nucleus, r = (Ze) 2 /(ajksT) is the conventional parameter that 
characterizes the strongness of the Coulomb interaction, and f(r) is given by 1 ^ 1 

/(r) rj 0.73317 - 0.39890r + 0.3414ir 1/4 + 0.05484r _1/4 . (C-19) 

The integration approximately gives for x = E u aj / (he) < 1 

IS- 1/(r) x 4 l 1 (/(r))2 x 6 1 (/(r))3 x 8 l 1 (/(r)) V ° 

(C-20) 
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-■ion/ 



To use this expression for the case of x > 1, we set the maximum value as {S- u 
max(l, (<Sion)) where (5; on ) = 1 corresponds to the case without the correction. 

C.2.3. Absorption on free neutrons 

The total cross section of the absorption of electron neutrinos on free neutrons 
is given 

X 1 W*\ (C-21) 



ai) / 1 + 3g\\ ( E v + AipV , / in r <- 



m P c^ 



1 



E u + A 



up 



where A np = m n c 2 — m p c 2 , and W^ h is the correction for weak magnetism and recoil 
of neutron. We use the exact expression derived by HorowitzP^HJ By contrast with 
the corrections in the scattering cross section on free nucleons, the corrections in 
the absorption do not show the bad behavior at low neutrino energies. Similarly, 
the total cross section of the absorption of electron anti-neutrinos on free protons is 
given by^ 1 

2 r / .2 



al) / 1 + 3g 2 A \ / E p - A np \ 1 ( m e c 
a = (7 ' 



l 



E A 



1 \ 4 J \ m e c A 

Again, we use the exact expression derived by Horowitz!^^ for . 



W; h . (C-22) 



p 
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